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ABSTRACT OF THE DISSERTATION 


Harmonic and Anharmonlc Properties of Diamond 
Structure Crystals with Application to the 
Calculation of the Thermal Expansion 
of Silicon 

by 

Keith Herbert Wanser 

Doctor of Philosophy in Physics 
University of California, Irvine, 1982 
Professor Richard F. Wallis, Chairman 

Silicon has interesting harmonic ar:d anharmonic pro- 
perties such as the low-lying transverse acoustic modes at 
the X and L points of the Brillouin zone, negative 
Griineisen parameters, negative thermal expansion and 
anomalous acoustic attenuation. In an attempt to under- 
stand these properties, a lattice dynamical model employing 
long-range, nonlocal, dipole-dipole interactions has been 
developed. Several interesting features of this intei ac- 
tion are found and discussed. Analytic expressions for the 
Griineisen parameters of several modes are presented. These 
expressions explain how the negative Griineisen parameters 
arise. It is found that short-range anharmonic ity is 


xiii 


» 


inadequate to explain the experimental data for the 
Gruneisen parameters. Application of this model to the 
calculation of the thermal expansion of silicon from 5K 
to 1700K is made. Good agreement with experiment is ob- 
tained from 17K to the melting point. The thermoelastic 
contribution to the acoustic attenuation of silicon is 
computed from 1-300K . Strong attenuation anomalies 
associated with negative thermal expansion ere found in the 
vicinity of 17K and 125K . Lxectrostatic multipole in- 
teraction energitJ ,n a nonlocal dielectric medium and the 
Ewald method and its application to dipole sums are dis- 
cussed in detail in the appendices. 


xiv 


INTRODUCTION 


The prevent study was motivated by a desire to under- 
stand the acoustic attenuation of silicon at low frequency 
and low temperature. Professor Joseph Weber reported 1 an 
interesting narrow dip in the mechanical quality factor or 
Q of a large, ( ~1 meter long) nearly perfect, single cry- 
stal of silicon in the vicinity of 12K at 3.4KHz. The 

2 

extremely large Q of such crystals, of the order of 

9 

10 , caused them to be of Interest as possible gravita- 

tional wave detectors or ultrastable clocks. 

. 3 

Existing anharmonlc calculations and measurements 
at 300-500 MHz show no peak in the acoustic attenuation 
of silicon at low temperature . Studies of impurity in- 

4 

duced. attenuat ion showed this mechanism to be unlikely, 
since the dip in Q was too narrow to be fit by a relaxa- 
tion type expression with an activated relaxation time. 

Furthermore, the boron impurity content of the silicon 

14 3 

crystal was 8x10 atoms/ cm , a rather low impurity con- 

centration. 

After considering several possible attenuation mechan- 
isms, the author became aware of the negative thermal ex- 
pansion of silicon. A calculation of the thermoelastic 
contribution to the acoustic attenuation of silicon showed 
strong attenuation anomalies near 13 K and 17K due to 
the change in sign of the thermal expansion coefficient. 
This calculation, together with earlier studies, suggests 


1 



that the anomaly in the Q is an anharmonic effect and 
not due to impurities. 

Since third order elasticity theory is unable to yield 

5 

the negative thermal expansion of silicon, it was realized 
that elasticity theory could give misleading results for 
anharmonic properties of diamond structure crystals. It 
was thus desireable to develop a realistic lattice dynamical 
model for harmonic and anharmonic properties of silicon. 

With this in mind, it was decided to calculate the thermal 
expansion of silicon. 

Silicon has negative Griineisen parameters for several 
modes and negative thermal expansion as a consequence of 

g 

this. The early work of Barron pointed out that an fee 
crystal with only nearest neighbor Hooke's Law central po- 
tential interactions would have a thermal expansion that is 
negative at all temperatures. Blackman' considered an 
ionic lattice model with the zinc blende structure and 
found negative Griineisen parameters in the elastic region. 
This work, while interesting, did not show how the negative 
mode gammas arise, and no analytic expressions for the mode 

g 

gammas where given. Bienenstock made a calculation of the 
thermal expansion of germanium but the approximations made 
were somewhat drastic and the physical mechanism for the 
negative mode gammas was not recognized. Also, no analytic 
expressions for the mode gammas were presented. Similar 

9 

criticisms apply to the calculations of Dolling and Cowley 


3 


and Jex, 10 aa discussed In detail in Chapter 2. 

In Chapter 1 we develop the harmonic lattice dynamical 
model, including the investigation of long-range nonlocal 
dipole interactions. In Chapter 2 we develop the enhar- 
monic model, which is a consistent extension of the har- 
monic model. Analytic expressions for the Gruneisen para- 
meters of several modes are presented and the model applied 
to a calculation of the thermal expansion of silicon. In 
Chapter 3 we present a calculation of the anomalous ther- 
moelastic attenuation in silicon and in Chapter 4 a sum- 
mary of conclusions is given. Several details pertaining 
to the discussion in Chapters 1 and 2 are given in the 
appendices . 
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Chapter 1 


Harmonic Lattice Dynamics of the 
Diamond Structure 


"0 Lord, how manifold are thy works! 
in wisdom hast thou made them all: 
the earth is full of fchy riches." 

Psalm 104:24 


INTRODUCTION 

In this chapter we develop a model for the harmonic 
lattice dynamics of diamond structure crystals. The chap- 
ter is divided into four main sections. First we discuss 
some elementary properties of the diamond structure. Se- 
cond, we develop a model for the short-range contributions 
to the dynamical matrix. Third, we develop a model for the 
long-range contribution to the dynamical matrix employing 
nonlocal dipoles. This is a lengthy section and requires 
several side discussions to be found in the appendices. 
Fourth, we present analytic expressions for the elastic 
constants and for phonon frequencies along symmetry direc- 
tions using the model developed in the previous two sec- 
tions. We also present fits to experimental data and 
phonon dispersion curves for silicon in this section. 


5 
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Figure 

tetrahedral 


1. The diamond crystal structure showing the 
bond arrangement. 
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Figure 2. The rhombohedral primitive cell of the 
face-centered cubic lattice. The primitive translation 
vectors *i’*2'*3 connect the lattice point at the origin 
with the lattice points at the face centers. Shown are 
conventions used in the present work. 



D IAMOND STRUCTURE 


The purpose of this section Is to define some conven- 
tions that will be used throughout this thesis, as well as 
acquaint the reader with some basic properties of the dia- 
mond structure. 

The diamond structure is a face-centered cubic space 
lattice with two atoms per unit cell.* (see Fig. 1) It may 
be visualized as two interpenetrating fee lattices with one 
shifted relative to the other by one quarter of the body 
diagonal of the conventional cube. The lattice position 
vectors are given by 

ft(iK) - ft(Z) + ft(K) , (1.1) 

with Bravals lattice vectors 

R(Z) - + i 2 a 2 + Z 3 a 3 , (1.2) 

and the basis vectors 

R(K) ■ ^(ii^+ag +a 3 ) for < - 0,1 . (1.3) 

Here ^1**2' *3 are arbitrar y integers. The primitive 
translation vectors a^ are given in terms of the carte- 
sian unit vectors by 

a l ' § ( ®1 + ®2 } d.4a) 

a 2 ■ "j( e 2 +e 3^ (1.4b) 

a 3-f ( W (1 * 4C) 


10 


11 


where a is the conventional cube edge which is 5 .43 ok 

2 

for silicon at 300K. (see Fig. 2) Thus we have in terms 
of cartesian coordinates 

it( XK ) - |[ + * 3 + K/2)e x + + i 2 + k/ 2)&2 

+ d 2 +i 3 +K/i')e 3 ] . (1.5) 

In the following we denote 

IU*|i'K'} - 3(fK) - fi(l'K') . (1.6) 


The volume of the primitive cell is 


n o ■ l‘V <a 2 xa 3 ) l ■ 


a 


(1.7) 


and the mass density of the perfect crystal is 


2M 8M 

p " a r 

o a 


( 1 . 8 ) 


where M is the mass of the atom. 

For the first few neighbors of a given atom we have 
the following. 

r - - - first neighbor distance 

o 4 


4 first neighbors. 


(1.9a) 


r 2 “ “ 1.6329932 r Q - second neighbor distance 


12 second neighbors. 


(1.9b) 


12 


r 3 • * ■ 1.9148542 r Q - third neighbor distance 

12 third neighbors, (1.9c) 

r^ ■ a - 2.3094011 r Q - fourth neighbor distance 

6 fourth neighbors. (1.9d) 

Thus we see that the fourth neighbor distance is the first 
to exceed twice the nearest neighbor distance. The tetra- 
hedral bond angle is given by (see section on angle bending) 

0 (O) - cos’* (-1/3) - 109.47122° . (1.10) 

This is the angle defined by a given atom and any two of 
its nearest neighbors, when in the equilibrium positions. 

We now discuss the reciprocal lattice. The reciprocal 
lattice to the face-centered cubic lattice is the body- 

3 

centered cubic lattice. The reciprocal lattice vectors 
are given by 

£ - mjSj + m 2 ^ 2 + n^bg (1.11) 

with m^m^ntj arbitrary integers. The primitive trans- 
lation vectors of the reciprocal lattice are given in 

terms of the cartesian unit vectors e^ by 


13 


b x - 'T^ 1 + ® 2 "®S ) ( 1 . 12 a) 

S 2 * T- ( -*l +S 2 +S 3 ) (1 - 12b) 

b 3 - ^(e x -e 2 +e 3 ) . ( 1 . 12 c) 

Thus wo have In terms of cartesian coordinates 

G - — [ (n^ - m 2 + “ 3 )®! + (“l + m 2 ” “3^2 

+ (-m^ + m 2 + nigjeg ] . (1.13) 

Just as we have a primitive cell in real space, we also 
have a primitive cell in reciprocal space. This minimum 
volume cell is called the first Brillouin zone. It is the 
Wigner-Seitz cell in reciprocal space. Figure 3 shows the 
first Brillouin zone of the face-centered cubic crystal 
lattice . 

At this point it is appropriate to make a brief com- 
ment on the stability of the diamond structure. Carbon, 
silicon, germanium and grey tin crystallize in the diamond 
structure. Th„ first-order semiconductor to metal transi- 
tion occurring in tin at 286K and atmospheric pressure is 
accompanied by a change in structure from the diamond cubic 

(grey, a-Sn) 0 the body-centered tetragonal (white, 

4 

0-Sn) form. This transition has been observed in silicon, 
germanium and a number of III-V compounds at high pres- 
sures. For silicon, a sudden drop in resistivity of over 
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14 


five orders of magnitude Is observed In the vicinity of 
120-150 kbar. 4 

Other crystal forms have also been found. A metasta- 
ble, semiconducting, hexagonal diamond ("Wurtzite" ) struc- 

g 

ture has been observed for carbon and silicon , as well as 

4 5 

a metallic, bod> -centered cubic structure for silicon. ' 

7 

Electronic structure calculations for silicon have shown 
that the diamond structure is the lowest in energy for a 
number of plausible crystal structures, (see Fig. 4) The 
next structure lowest in energy is the hexagonal diamond 
and then the P-tin structure. By way of contrast, above 
~-100°C ice Ih , which is the hexagonal diamond struc- 
ture for the oxygens, is more stable than ice Ic , which 

8 

is the ordinary diamond structure. It appears that the 
two structures for ice are very close in energy, as they 
are for silicon as shcv'n in Fig. 4. The reason for this 
is that the tetrahedral coordination is preserved and the 
nearest neighbor distance is essentially the same for the 
two structures. The second neighbor distances are also 
the same, but the orientation of the second neighbor posi- 
tion vectors differ between the two structures. 
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cubic 


Figure 3. First Brillouin zone of the face-centered 
lattice with symmetry points and lines labeled. 
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Figure 4. The diamond, hexagonal diamond, 0-tin , 
hep , bee , and fee structural energies (in units of 
Ry/atom) as a function of the atomic volume (normalized to 
the measured free volume) for Si . The dash d line is the 
common tangent of the energy curves for the diamond and the 
0-tin structures. (Courtesy M.T. Yin, Ref. 7) 






HARMONIC APPROXIMATION AND EQUATIONS OF MOTION 


Before constructing the harmonic model, a few words 
concerning the harmonic approximation and equations of 
motion are in order. The harmonic hamiltonian is 

H, - T + « 2 (1.14) 

where T is the kinetic energy and $ 2 is the contribu- 
tion to the crystal potential energy which is quadratic in 
the nuclear displacements. Thus we have for the diamond 
structure 

P^(jtK) 

T “ ,? SM ' (1.15) 

IKOL 

and 

*2*1 S )U 0 (i 'K ' ) , (1.16) 

1V( f 


with 


$ a 0 (£KjrK') 


jTi 

du (iK) du a n 'X v ) 

d p 


(1.17) 


Here we have denoted P (i-K) and u (IK) as the a. 

a a 

cartesian component of the momentum and displacement, re- 
spectively, of atom (£K) . Note that u^(Zk) is mea- 
sured from the rest position specified by Eq. (1.1). The 
* p (ZK 1 1 'k ') are known as the second order atomic force 
constants. 4 is the potential energy of the crystal as a 
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function of the nuclear displacements and the derivatives 
in Eq. (1.17) are to be evaluated with all atoms in the 
rest positions. Using the Hamiltonian Eq. (1.14), one can 
show that the equations of motion are 


Mu (IK) 
CL 


« ae (XK|i'K')u g U'K') . 


(1.18) 


This equation holds classically and also quantum mechani- 
cally, as can be verified by using the equal time commuta- 
tion relations for P^dK) and u^Uk) as well as the 
Heisenberg equation of motion. 

If we now seek normal mode solutions of the form 


where k 
Eq. (1.18) 


W (K ) , 7 * j* . 

u <«) - 

o. JU * 

is the waveveotor and u) is the frequency 
becomes 


(1.19) 
, then 


a) 2 W (K) - Z) D „ (KK ' |£)W S (K ') , (1.20) 

CL |3K' 0 

where D a g(KK'|k) is the dynamical matrix defined by 

D fl (KK'jjc) • i £ « (iK |i'K Oe - ^' (1.21) 

Note that in this thesis we will be using the first form 

of the dynamical matrix, Eq. (1.21), as opposed to the 

g 

second form of the dynamical matrix. The two differ by a 
k dependent phase factor. 
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Using the invariance of the force constants under a 
crystal lattice translation it(m) 10 we have 

* g (i +m,K \l' + m,K ') - (IK \ l*K* ) . (1.22) 

Letting m - -I in Eq. (1.22) gives 


* afi (0,K\t' - « ap (XK|i'K') . 


(1.23) 


Thus Eq. (1.21) can be written as 


D aS (K«'|k) 


■ ill (0,K|i'-/,K')e li ‘ ,[?C1 ' ) " R(1)1 . 
a * 


(1.24) 


Neglecting surface effects, we can relabel the sum in Eq. 

(1.24) so that 


D a e< KI< ' 


k) 


i T, 

M „ 
m 


‘ a p (c ’ K 


m, K ' )e 


ik*R (m) 


(1.25) 


From Eq. (1.25) we see that the dynamical matrix is inde- 
pendent of the unit cell index l in Eq. (1.21). 

Since the order of partial differentiation is inter- 
changeable, we have from Eq. (1.17) the "flip” symmetry 
property 


• p (jUC|i'iC') - | JtK) . (1.26) 


Using this fact plus translational symmetry, Eq. (1.22), 
one has the Hermit ian property 
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D ag (KK'|?) - D* a (K'K|ic) , (1.27) 

where the * denotes complex conjugate. From the defini- 
tion Eq. (1.21) we also have 

D*g(Ktc'|is) - D ap OCK'! -it) . (1.28) 

Now for the diamond structure, D _ ( k K ' | k ) is a 
6x8 Hermitian matrix. Thus it will have six real eigen- 
values for each value of it . The eigenvalues we label 
2 -* 

w (kj) , where j - 1 , 2 , • • • ,6 denotes the branch index. 
Furthermore, to each value of it we can find a complete 
orthonormal set of eigenvectors ( k { lc J ) that obey the 
eigenvalue equation 

oj 2 (kj)e a (K |itj) - D aj3 (KK'|it)e 3 (K'|itj) , (1.29) 

0 K 

the orthonormality property 

£ e * ( k | it j ) e ( K | it j ' ) - 6 , , (1.30) 

and the completeness relation 

L e a (K|itj)e g (K'|kj) » 6^6^, . (1.31) 

J 

Without loss of generality we also have the following pro- 
perties, which correspond to a particular set of phase and 

1 9 

branch labeling conventions 
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aj(icj) ■ , (1.32) 

e o (Kj£j) - e * ( K | -k j ) . (1*33) 

With these preliminaries out of the way, we can de- 
rive some properties of the force constants and dynamical 
matrix imposed by inversion symmetry plus translation in- 
variance . 

Now consider the effect of a crystal symmetry opera- 
tion represented by (S|7(S) +#(m)} . Applied to the posi- 
tion vector given by Eq. (1.1), this operation transforms 
it as 

(S | v (S) +jt(m)}ftUK) - 3»R(JU) + v (S) + R(m) 

s it (LK) . (1.34) 

S is a 3x3 real orthogonal matrix representation of 

one of the proper or improper rotations of the point group 

of the space group, v(S) is a vector smaller than any 

primitive translation vector of the crystal, and it(ra) is 

one of the translation vectors, Eq. (1.2). For a review 

of crystal symmetry operations, we refer the reader to the 

13 14 

article by Maradudin and Vosko and the book by Lax. 
Furthermore, the invariance of the potential energy under 
the operation {S|v(S) +Ff(m)} produces the transformation 
law for the second order force constants 


I 
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«^<LK|I/K') - Z S^l^dKirK') . 
For the inversion operation we have*'* 

S a0 " “ 6 a0 

where 6^ is the Kronecker delta symbol, 
v(S) - R (1 ) - a/4(5 1 +e 2 +e 3 ) , 


and 


5(m) - 0 . 


(1.35) 


(1.36a) 


(1.36b) 


(1,36c) 


Thus from Eq. (1.34) the position vector is taken into 

it (LX) = -5(2 K) + 5(1 ) (1.37) 

so that in particular 

5(L,1) - -5(2,0) + 5(1) - 5(-2,l) (1.38a) 

S(L,0) - -5(2,1) + 5(1) - 5(-2,0) . (1.38b) 

Therefore the effect of the inversion operation is to in- 
terchange the sublattice labels K -> 0,1 and take 5(2) 
into - 5(2) . This point will be important later in the 
section on nonlocal dipoles. Thus for the inversion opera- 
tion applied to the force constants, Eq. (1.3.5) becomes 


* ae aK| L 'K') - 

t ag UK|i'K') , 

(1.39) 

or more explicitly using Eq. 

(1.38) 


«afj<*.0K,0> - 

» ae (-l,l|-l'.l) 

(1.40a) 

» ag (!,0U',l) - 

‘^(-1,1 1-1', °) 

(1.40b) 


Schematically, we can write the dynamical matrix. 


Eq. (1.21), in block form 


[d^Ock'IS)] - 


D ae (°,0|?) 

D^TTTffTTT 


D (0, 1 (k) 


«3^' v 

The Hermit ian property, Eq. (1.27), is 




D a e<°.°l i ‘> - DB o «>.0|f> 


D oe a,i|S) 

D ae (0,l|J) 


Consider 


D po (l,l|k) 

Dj a <l.<>|S) • 


D^a.iiS) 


- I P ^U.lU'.De 


-ilc* [ftd)-ft(X') j 


using Eq. (1.40d) we have 

D „B -S p 

JL 


using lattice translation symmetry, Eq. (1.22), and 
ing the summation variables we obtain 

o oS a.i|S) - D* s (o,o|iT) . 

Consider next 

l 


(1.41) 

(1.42a) 

(1.42b) 

(1.42c) 

(1.43a) 

)] 

(1.43b) 

relabel- 

(1.44) 

] f 

(1.45) 
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t 


us lng Eq . ( 1 . 40c ) we have 


D aS a,°|E) - i 2 • ae (-l,0|-l'.l)e- lJ,[ * (1, - 5(i ' )1 


(1.46) 


again using lattice translation symmetry plus relabeling 
the summation variables we obtain 

D a3 (l,0|£) - D* g (0,l|k) . (1.47) 

Combining Eqs. (1.42c) and (1.47) we further obtain that 


D a9 (0,l|k) - D po (0,l|k) . 


(1.48) 

In summary, the form of the dynamical matrix is reduced to 

, (1.49) 


(D aj 3 (KK'| k )] - 


D a 0 (O,O|£) 

D aB (°,l|E) 

D^(°,l|k) 



with 5(0, Oj?) Hermitian and D(0,l|it) symmetric." 6 This 


form holds for any lattice with two like atoms per unit 
cell that are interchanged by the inversion operation. 


CONTRIBUTIONS TO THE DYNAMICAL MATRIX 


Since the potential energy can be written aa a sum of 
contributions trcj various types of interactions, the force 
constants Eq. (1.17) and the dynamical matrix Eq. (1.21) 
can also be written as a sum of contributions from the 
various interactions. This allows separate computation of 
the contribution to the dynamical matrix from each type of 
interaction, and the total dynamical matrix can then be 
obtained by summing the different contributions. 

For the model considered here, we divide the potential 
energy into a short-range part and a long-range dipole- 
dipole part so that 


* - 4 s - r - + S dd 
2 2 2 * 


(1.50) 


and we model the short-range contribution by 


.S.R. _ .1° . s 2° , .3° . -4° .(0) 

*2 *2 + *2 + S 2 + *2 + *2 


(1.51) 


1° 2° 3° 4° 

Here 4 2 , , * 2 and * 2 are the harmonic contribu- 

tions to the potential energy from first, second, third 

(6 ) 

and fourth neighbor central potential interactions, $ 2 
is the harmonic contribution from nearest neighbor angle 
bending and $ 2 is the harmonic part of the long-range, 
nonlocal dipole- dipole interaction energy. 

Using the same notation, we label the contributions 
to the dynamical matrix 
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D ae (KK '|£> - D^; R (kic'|£) + D^Ocic'lfc) , (1.52) 


with 


0®j*-(KK'|S) - D^OdC'lf) +Df„(KK'|S) + d£<KK'|S) 


ad 


QL& 


«e 


+ D^(KK'fk) + } (K K ' |k) . 


atfi 


(1.53) 


We now proceed to a discussion of the various terms in Eq 
(1.53) and then the dipole-dipole part of Eq. (1.52). 


CENTRAL POTENTIAL CONTRIBUTIONS TO THE DYNAMICAL MATRIX 


For a two body central potential Interaction, the po- 
tential only depends on the magnitude of the separation of 
the atoms. Denoting the potential for an atom of type k 
interacting with an atom of typo tc ' at a distance r 
from it by ® KK# (r) , we can write the crystal potential 
energy due to central potentials as 

•-is £' 0 ,(|R(1K|X 'K') +u(iK |i'K') |) , (1.54) 

C 2 IK t'K' 

where 

U(ZKji'K') - u(iK) - u(i'K') . (1.55) 


The factor of 1/2 in front of the sum compensates for 
the fact that each interaction is counted twice in the sum, 
the prime on the second sum means exclude terms with 
(IK) - U'K') . 

Expanding the potential (r) in powers of the dis- 

ft ft 

placements leads to 


a KK ,()3(XK \l'K') + u(£K | f ' 1C ' ) |) - 0^ ,(|3(2K |rO |) 


K K 


+ Ej (IK \t 'k')u (iK \l'K') 

Of a ® 

+ | £ ® ag (XK [X'lC')U a (XK |£'K')u g (iK Ji'K') 


+ | £ ® a ^ Y UK|i'K')U a (iK|jl'K')u p (XK|f'K')U Y (iK| ( e'K') 

aP y 


+ (higher order in displacements). 


(1.56) 
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with the expansion coefficients given by 


17 




— 2 ® * , (x* ) 
r w kk v ' 


r-R (IK 1 1 'k ') 


i X x 0 

e^dKH'K') - j-Sji 


® #(r) ® * , (r ) 

*KK VA ' r W KK ' Vi ' 


+ »KK' (r > 


?-£(flC I'K' 


afly 


UkU'O - ! ^£g 


1 r 3 


3®" , (r ) 3®' , 

,_x ®KK vi 7 , ° V KK 

®KK' (r) T + 2“ 


x 6 -hx_6 +x 6 ,, 

SL P_av y fl g 


®^' <r) *F ®KK' (r> 


r-£(XK £ 


and 


? - 2 «„#„ • 
a 


Here we have denoted the derivatives by 


<t> ' (r) 

®"(r) 

/// , x 

® (r) 


d® (r ) 
dr 

d ® ^r) 


dr 

£^0rl 

dr 3 


Using Eqs . (1.54) and (1.56) we obtain the harmonic 


(1.57a) 


) (1.57b) 
(r) I 


'k ') 
(1.57c) 

(1 . 57d ) 

(1.58a) 

(1.58b) 

(1.58c) 

contri- 


bution to the crystal potential from central potentials 
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i 


2c 


i £ £ 0 fl (XK|i'K')u (iK|Z'K')u a UK X'K') . 

4 IK a i'tc'9 ap ° P 

(1.59) 


We now proceed to compute this for first, second, third 
and fourth neighbors. 


FIRST NEIGHBOR CONTRIBUTION 
Writing 

1 ° 2 ° 2 ° 4 ° 

* 2 c “ *2 + *2 + *2 + *2 (1 ^ 0) 

we have upon evaluating the sum in Eq. (1.59) over nearest 
neighbors 

1° 4 r 

*2 -i s 2 0 0 (X,O|X+6 l)u(X,0|X+6 l)u„(X,0|X+6.,l) 

45 Xotf i**l(_ 1 “ i P i 

+ 0 ag (X, 1 |X+ 6 i + 4 , o ) u a (X, 1 |X+ 6 i+ 4 ,O) u p (X,l|X+ 6 i+ 4 ,O) , 

(1.61) 

where 

“ 0 > ^2 "” a j. ’ ^3 * “ a 2 * "^4 = " a 3 (1.62a) 

and 

t i+4 - for 1 < i < 4 . (1.62b) 

The are the translation vectors to the unit cells of 

the nearest neighbors from the atom of interest. From the 
definition Eq. (1.57b) we have the lattice translation 
property 

® ag (X + m,K|X' +m,ic') - (-* K |X'iO , (1.63) 

and the further property that 

(Xk|X'k') * ® a p(X'K'|4K) . (1.64) 
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Thus we have 


33 


» ap (l f O|l +a i# l) - 0 ap (O,O| 6 i ,l) (1.65a) 

® a$ (l,l|X + 6 i+ 4 ,0) - « ap (0,l|6 1 + 4 ,0) , (1.65b) 


and further 


® ap (l,l|l+6 i+4i 0) - 0 ag (O,O|6 i ,l) . (1.66) 


Thus we can write Eq. (1.61) as 


2° -* A S »ae <0 >°l 6 i- 1) 


fa3 i-1 


u a (i,O|X+ 6 i ,l)u g U, 0 |i+ 6 i ,l) 


+u (i,l|l-6 i ,0)u p (l,l|l-6 if 0) 


(1.67) 


From the definition Eq . (1.57b) we can directly compute the 
four matrices (0, 0 | 6 1) . We exhibit them explicitly 


0 (0, 0 | , 1 ) 


0 (0,0|6 3 ,1) 


/ a 

3 

3 \ 

/ a 

3 

-3 \ 

0 

a 

e 

ST(0,0|6- f l) - 3 

a 

-3 

\ 3 

3 

a / 

\-3 

-3 

a / 


/a -6 -3 \ / a -3 M 

-3 a 3 ®(o,o|6 l) - -3 a -3 

\ -3 3 a / \ 3 -3 a / 


( 1 . 68 ) 


where we have defined the force constants a and 


a 


I ®i (r o ) + sf; ®i (r o ) 

5 »i (r o ) - 57- ®i <r 0 ) 

O 


3 as 


(1.69a) 

(1.69b) 
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For the sake of being explicit, we have labeled the near- 
est neighbor potential function $j.(r) . The derivatives 
are to be evaluated at the nearest neighbor distance Eq. 
(1.9a). 

Since Eq. (1.18) can also be written as 


Mli UK) 
a 


a *2 

9u^?Tk7 ' 


(1.70) 


if we now substitute the normal Tiode solution Eq , (1.19) 
into Eq . (1.20) we obtain the relation 


E D 

0K ' 




(kk' |ic)w g (O 



e ~ik' RU )+iurt 


(1.71) 


where it is understood that we have substituted the normal 
mode solution after performing the derivative. Direct com- 
putation using Eq. (1.67) yields 

3*2° 4 

au g TxT0) ~ ^^(O'Ol 6 !’! )u p(^0|^6i,l) (1.72a) 


di 2 


E E « (0,0(6 ,l)u (1,1 

i-1 13 “ p 1 p 


1 - 6 ^ 0 ) 


(1.72b) 


Using Eqs . (1.71) and (1.72) we obtain 


d£<0.0|S) E t 0(0,015.,!) , (1.73a) 


i=l 
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D*°(0,1|£) - - -| ^ og (0,0j5 i ,l)e 1 ^ ^ . (1.73b) 


Using the matrices Eq. (1.68), we have more explicitly 




•♦1 0 . 

d a (0,0|k) 


0 

1 

0 


(1.73c) 


SECOND NEIGHBOR CONTRIBUTION 


Evaluating the sum in Eq. (1.59) over second nearest 
neighbors we obtain 

0 a0 (X,O|X+u i ,O)u a (i,O|i-Hi i ,O)u g U,O|l+u i ,O) 
+ ^ +tJ ’i+12 , ^^ u a^ , ^‘ ^ , * + ^i+12 , ^^ u p ^ +M, i+12 , '*‘^ * 

■m 

(1.74) 

where 

•> -► -> 

14 1 “ a l » ^2 “ a 2 ’ ^3 ’ a 3 ’ 

-» .-+ •* . .-*• -*• . -> . 

^4 " a 3~ a 2^ ’ ^5 “ ^ a l“ a 3^ * ^6 ” ^ a l” a 2 ' • 

M- 7 “ ” M 1 , M-g ” -p- 2 * M-g " -M- 3 , (1.75a) 

^10 " “^4 ’ 4 11 " ~^5 * ^12 " ~^6 * 

and 

^i+12 * for 1 < 1 < 12 • (1.75b) 

The ^ are the position vectors to the second near- 
est neighbors from the atom of interest. In this case, 
the for 1 < i < 12 are the position vectors from 

the origin to the twelve nearest face centered sites. 

Using Eqs. (1.63), (1.64) and (1.75b) we have 

(i, o U+Ui» 0) =* ® a g (0, 0(^,0) , (1.76a) 


. 12 
-itS 

* ictf i-1 
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» U.iU-u^.u - « .(0,01^, 0) 


a 


So we can write Eq. (1.74) as 


12 


•2 -kL ^Ads' 0 ’ 01 ^’^ 


la3 i-1 


u a (i,0|l+^ i ,0)UgU,0|l+u 1 


+ u a (X,l|l-u i ,l)u g U,l|l-^ i ,l) 


(1 


Noting that 


i+6 


, 1 < i < 6 


(1 


we obtain 


® a g (0, °^i+6' 0) " 0 tt p (0 » 0 1^1,°) , 1 < i < 6 . (1 

Therefore we only need compute six matrices. Using Eq 
(1.57b) one obtains 


$ (0,0 |^ 1 , 0) 


O' 
0 

0 0 X , 


'u v 

V n 


®(o,o|n q ,o) - /u o 

oxo 
\V o v.) 


0 (0, 0 I Hg, o) 


$ (0, 0 I m > 4 , 0) 


/X 0 0 

( 0 u V 

\ 0 V n 

/ [i -v O' 

\-V M. 0 

\ 0 0 \> 


0 (0, 0 lug, 0) 


/X 0 O' 
I 0 M, -V 
\ 0 —u 


0 (0, 0 In-, 0) - 


( \L 0 -1^ 

oxo 

\-V 0 \l) 


(1 


Here we have defined the force constants and X 


.76b) 

, 0 ) 

.77) 

.78) 

.79) 


80) 

as 



t 
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u 


V 


l(^2 (r 2 ) + 
l(®2 (r 2 ) " 


^2 (r 2 } 
r 2 

• a<r 2 ) 


X 


(u-iO 


®2^ r 2^ 


(1.81a) 


(1.81b) 


(1.81c) 


Again, for the sake of being explicit, we have labeled the 
second neighbor potential function The derivatives 

are to be evaluated at the second nearest neighbor distance 
Eq. (1.9b). Direct computation using Eq. (1.77) yields 
2 ° 

12 

aT777oT- ESv (, ' , K- 0)u )“' 0| “*i' 0) (1 ' 82a) 

Of 1™1 P 

2 ° 

)ii 12 

a TT O TTl ~ 1 ^'«p ( 0 '°K l ”V 1 ’ 1 , 1 V ) ' (1 ' 82b) 

Using Eqs . (1.71) and (1.32) we obtain 



(0,1 fk) 


0 



( 0,0 ( it ) 



6 

s < 

i-1 


aP 


( 0, 0 |p. ± , 0 ) r 1 -costMTj)] 


(1.83a) 

(1.83b) 


The first equation expresses the fact that the second 
neighbor interaction produces no coupling between the 
< * 0 and k m 1 fee sublattices. In obtaining the second 
equation we have also used Eqs. (1.78) and (1.79). 



THIRD NEIGHBOR CONTRIBUTION 


Evaluating the sum in Eq. (1.59) over third nearest 
neighbors we obtain 

* ag U,0|x+r 1> l)u a (i,0|X-»-T 1# l)u e U,0|l+T 1 ,l) 

+ ®a£ ^ i+12* ® ^ u ot ^ » * I ^ +T i+12' ^ I i+12' ^ 

(1.84) 

where 


T 1 

. -* . -+■ 
( a 3 ” a l ) * T 2 

(a 3 - a- 1 - a 2 ) , t 3 - (a 3 - a g ) 

9 

-► 

T 4 

• (a 2- a l ) * T 5 

^ a 2 ” a 3 ' T 6 ’ ^ a 2 ■ a i “ a 3 ) 

9 

-+• 

T 7 

^ a l ” a 2 " a 3 ^ 

’ t 8 “ ^ a l ” a 2^ * 7 9 ^ a l “ a 3 ^ 

9 

T 10 

* "< a 2 + a 3 } '"ll 

.-*■ -+ . -* -► . 
- -(a x +a 3 ) , t 12 - -(a x + a 2 ) 

9 




(1.85a) 


T i+12 

-t for 1 < i < 12 . 

(1.85b) 


The T i are the translation vectors to the unit cells of 
the third nearest neighbors from the atom of interest. 
Using Eqs . (1.63), (1.64) and (1.85b) we have 


^ > 0 | i-HT ^ , 1 ) » 

i^.i) 

(1.86a) 

(i » 1 i i+r i+12 ,0) * ® a g (0 '°J 

iTi-D • 

(1.86b) 


i 12 

i £ 2 

* tat 1-1 


39 


M- 


40 



5' n' V 


i' U' n' 


$(0,0 T 1t ,l) - /\' b' 6 ' \ $(0,0 T..,l) - /»' 6' »/' 


6' X' 6 

V’ 6' M, 



( 1 . 88 ) 

defined the force constants 

as 

n: [ ®3 (r 3 ) »3 <r 3 >1 

(1.89a) 

TT (, 3 (r 3 ) ®3 (r 3 )| 

(1.89b) 

3^ - Yj[»3<r 3 ) #3 (r 3 ) ] 

(1 .89c) 

(* +* u ) " [ TT ®3 (r 3 ) + llr 0 

® 3 ( r 3 ) ] • ( 1 . 8 9d ) 


Again, to be explicit, we have labeled the third neighbor 
potential function $^(r) . The derivatives are to be 
evaluated at the third neighbor distance Eq . (1.9c). 
Direct computation using Eq. (1.87) yields 


£ 2 $ ap (0,0|T 1 ,l)u ( J e,0|l+T 1 ,l) (1.90a) 

i“l 0 


2 

rn7 


£ S $ ae (0,0)T i ,l)u g (X,l!X-T i ,0) (1.90b) 

i”l 3 


Using Eqs . (1.71) and (1.90) we obtain 


D fp ( 0 -°I S > ‘ 3 £ ®aS ( 0 '°l T t> 1) 


(1.91a) 
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D ae (0 ’ 1 '^ ) * 


i 12 lk»T. 

* M 1?! ®a9 <0 ’ 0|T i’ 1)e • (L91b> 


Addition of the matrices Eq . (1.88) ^ives more explicitly 


D 3 ° (0, 0 | Ic ) - 1 o°) . 


(1.91c) 


FOURTH NEIGHBOR CONTRIBUTION 

Evaluating the sum In Eq. (1.59) over fourth nearest 
neighbors we obtain 



L L 0 (X, 0 |i+\ ,0)u (i,0|l+X i ,0)u B (jt,0|i+X i ,0) 

* lag 1-1 ap 

m 

+ ^ ^ ^ i+6 » ^ ^ I i+6 # ^ ^ I i+6^ ^ ^ 


(1.92) 


where 



ae 2 ' ^3 " ae 3 


-► -► -► 
-x 2 i ” "^3 

(1.93a) 

1 < i < 6 . 

(1.93b) 


The \ i are the position vectors to the fourth nearest 
neighbors from the atom of interest. In this case, the 
\ i for 1 < i < 6 are the position vectors from the 
origin to the six nearest simple cubic lattice sites. 

Noting that 

0 ap (l,O|i+\ i ,O) = ^ (0,0|Xi, ) (1.94a) 

$ aB (i, l|i+\ i+6 ,l) - ® a p(0,Ol\ i ,0) (1.94b) 

we can write Eq. (1.92) as 
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»f v<2 S 4 ae <0,0 ^i ,0) 


ictf i-1 


u a (i,0|i+X iJ 0)u g (l,0|X+X i ,0) 


+ u 0£ (i,l|i-X 1 ,l)u p (Z,l|i-X i ,l) 


(1.95) 


Since from Eqs . (1.93a), (1.63) and (1.64) 

<e ae ( °*°l X i+3 ,0) " <0 »°l x i' 0) » 1 £ 1 1 3 (1-96) 

we only need compute three matrices. Using Eq . (1.57b) 
one obtains 


© (0, 0 j x x , 0) - I\ u 0 OX £(0,0|X 2 ,0) - /\i u 0 0 


0 n 0 


0 X 0 


iO 0 p, ' 


^0 0 p, 


®( 0 , 0 (x 3 , 0 ) - 0 0 


0 a 0 


k 0 0 x 


(1.97) 


Here we have defined the force constants p. * and X" as 


«>4 (a) 


(1.98a) 


// // 


- 04 (a) . 


(1.98b) 


The fourth neighbor potential function we have labeled 
® 4 (r) and the derivatives are to be evaluated at the 
fourth neighbor distance, which is the conventional cube 
edge, (see Fig. 2) 
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One comment is appropriate at this point. Our 

tion for the force constant*, Eqs. (1.69, 1.81, 1.89, 

18 

has been chosen to agree with the paper by Herman. 
Direct computation using Eq. (1.95) yields 


aTTZTffy 


- £ £ 0 fl (0,0|\ ,0)u fl (1.0|l+\ ,0) 

i-l 0 aP ip i 


3*2 

3u a U, 1) 


6 

£ ® ag (o,o|\ i ,o)u g u,i!i-\ i ,i) 


Using Eqs. (1.71) and (1.98) we obtain 


D’ B < 0 ,l.|k) - 0 


2 

D ^(°,°|ff) - | 2 « aB (0.0|X 1 ,0)(l-coBg.? t ) 


Letting 


k - k 1 e l + k 2 e 2 + k 3 e 3 , 


we obtain the more explicit form 


D a°9 (0 ' 0|S) ‘ D al (0 ’°l S)6 a e ' 


with 


D’i(°,0|k) - | 




nota- 

1.98), 

(1.99a) 

(1.99b) 

(1.100a) 

(1.100b) 

( 1 . 101 ) 

( 1 . 102 ) 

)] 

(1.103a) 
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D^CO.Ojit) 

4° 

D33 ( 0 , 0 ( k ) 


,' S ln 2 (^) + r. in 2 (^)^- sln 2(V)' 


(1.103b) 


« . 2 
\i sin 


(¥) 


v . 2 k 2 a . » . 2 

-Ha sin ^ + ^- sin 


(¥)' 


(1.103b) 


Eq. (1.100a) expresses the fact that the fourth neighbor 
interaction produces no coupling between the k - 0 and 
K - 1 sublattices. In obtaining Eq. (1.100b) we have also 
used Eqs . fl.93a) and (1.96). 



ANGLE BENDING CONTE I BUT ION 


The angle bending contribution to the potential energy 
is a three body force. The reason for this is that three 
points are necessary to define an angle. Because of this 
feature, the angle bending contribution to the dynamical 
matrix requires much more computational effort to obtain 
than the central potential contributions. For reference 
purposes we include a more complete discussion of angle 
bending in Appendix A. In this section we only outline the 
results found there. 

We consider "pure" angle bending of the form 

Ae 2 (z,o|x+6 1 ,i|i+5 2 ,i)+Ae 2 (i,o|i+& 1 ,i|i+& 3 ,i) 
+A0 2 (l,O|i+6 1 ,l |X+6 4 ,1)+A0 2 U,O| f+6 2 ,l |i+6 1) 

+A0 2 (X,O| Z+6 2 ,l J X + 6 4 ,1)+A0 2 U,O|£ + 6 3 ,1| £+6^1) 

+A 0 2 ( £ , 1 |£+6 5 ,O|i+6 6 ,O)+A0 2 U,l | i+6 5 , 0 | i+6 , , 0 ) 
+A0 2 U,1[£+6 5 ,O|£+6 8 ,O)+A0 2 U,1| £+6 g , 0 | 1+6 ? , 0 ) 
+A0 2 U,1 j£+6 6 ,O|£+6 8 ,O)+A0 2 (£,l | £+6 ? , 0 | £ + 6 g , 0) 

(1.104) 

The are given by Eqs . (1.62). Here we have defined 

A0(£K | i. k | £ / k * ) as the change in angle between the atoms 
(i<) , (£' <■') and ( ) , having (£k) at the vertex, 
produced by the respective lattice displacements, (see Fig. 
5) The form of Eq. (1.104) is appropriate to "nearest 


ar 

8 ^£ 
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Figure 5. The angle 8 ( 4K | 4 * K * j 4 " k ” ) between any 
three atoms in the crystal. Atom (4k) is at the vertex 


of the angle. 
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neighbor" angle bending. There are six tetrahedral angles 
per a:om and two atoms per unit cell, thus twelve terms per 
unit cell. 

At this point we would like to make a few comments. 

There is some confusion in the literature concerning what 

is meant by nearest neighbor noncentral force. As we will 

see, our "nearest neighbor" angle bending potential will 

couple atoms that are first neighbors as well as atoms that 

are second neighbors, so in this sense it is not a purely 

nearest neighbor interaction. This point was noted by 
19 

Herman. However there is a difference between central 
potentials and what are termed central forces. We may see 
this by considering the force on an atom of type k =* 0 
produced by a variational displacement of its nearest 
neighbor in the same unit cell, all other displacements 
being zero. We have from Eq. (1.72a) 

1 ° 

1 ° d4 2 

F a (1 - 0) -~ au a U,o) -Z« ag (0,0|0.1)u B (i!,l) . (1.105) 
The line joining the center of the two atoms is 

£(1) "f(e 1 +e 2 + ®3 ) ' (1.106) 

and a displacement parallel to the line of centers is 

A 

u(X,l) - u n (e x +e 2 +e 3 ) . (1.107) 


This displacement produces the force 
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f 1 * (1,0) - (a + 20 )u |( (e 1 +e 2 +e 3 ) (1.108a) 

or 

f 1 * (l, 0) - 0^(r Q )u „ (e x +e 2 +e 3 ) , (1.108b) 

where we have used Eqs . (1.68) and (1.69). 

A displacement perpendicular to the line of centers 
Is 

u(£,l) - u x (e 2 -e 3 ) , (1.109) 

as can be easily seen since 

u x (e 2 - e 3 )»3(l) - 0 . (1.110) 

For this displacement we obtain 

f lO U,0) - (a-0)u x (e 2 -e 3 ) , (1.111a) 

or 

-►1° < r n > 

F U, 0 ) - v - u,(e 9 -e,) , (1.111b) 

r 0 1 4 J 

where again we have used Eqs. (1.68) and (1.69). For a 

"central” force, a displacement perpendicular to the line 

of centers of two atoms must produce no force. For the 

central potential to be a "central" force, we see from 

Eq. (1.111a) that this requires a ® 0 . This agrees with 

20 

the form of the matrices Eq. (1.68) exhibited by Herman 
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for the case of central forces. If 0{(r Q ) ^ 0 we see 
from Eq. (1.111b) that the central potential can produce a 
noncentral force. If only nearest neighbor central po- 
tential interactions were present, the condition of mini- 
mum potential energy would require 0^(r o ) “ 0 • kut in 
this case the lattice is unstable, as we will see later. 

If in addition second neighbor central potential interac- 
tions are present, the condition of minimum potential 
energy does not require 0^(r Q ) ~ 0 * Provided Q ) is 

not zero, we have an example of a first neighbor nonceutral 

force. The somewhat misleading statement has been made by 
21 

Keating that, ’’there are no noncentral purely first- 

neighbor interactions present in any nonmetallic crystal.” 

22 

Concerning this Ludwig has commented, "thif; is obviously 
not correct.’’ We also note at this point that the non- 
central force used by Keating is not "pure" angle bending. 
His noncentral coefficient 3 enters into the expression 
for the bulk modulus, whereas "pure" angle bending does 
not contribute to the bulk modulus. This is so because an 
isotropic homogeneous deformation leaves all angles in the 
crystal unchanged, and therefore any potential energy which 
only involves changes in angles will be zero under this 
type of deformation. The potential used here, Eq. (1.104), 
is an example of "pure" angle bending. 

Now to be more precise 
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Ae Uk|x'k' |x'k' 


e(iKji'K'fl'K # ) - 0 (XK I I'K ' | 1*k‘) , 


( 1 . 112 ) 


where 


COB 0 (IK | X'tC ' | i*K*) 


X(1K U'k')»x(AK I A*K*) 

{x(1k |A'k')||x(1k| J*K*) | 

(1.113) 


and 

x(XK)X'K') - XUK) - (1.114a) 

x(XK) - 3(Xk) + u(iK) . (1.114b) 

Here |x| denotes the magnitude of the vector x . The 
equilibi’ium angles are given by Eq. (1.113) when all atoms 
are in the rest positions 

COS 8 (0, (iK|i'K'M'K') - *i**l*'*’)'*U*l*'*‘) . 

Jlt(XK |X # K ') j [lt(XK \l"K") | 

(1.115) 


Note the symmetry property which follows from Eq. (1.113) 

cos 6 (iK |i'tc' | z"k") - cos 0(XK \1"k" \1'k') . (1.116) 

From Eqs . (1.112-1.115) it is clear that in general 
A6 (XK | i '< ' I l u K " ) is a nonlinear function of the displace- 
ments u^(Xk) . Thus our potential Eq. (1.104) contains 
anaharmonic terms. Provided that Ju(Xk|x'k'| -r [H(XK |x'k") | 


is small compared to one, we can obtain an expression for 


the change in angle, correct to first order in the dis- 
placements 
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A9(1K i'A'k'I I'k") - £ T) (ZK |Z'k' | i'*')u (lK|l'K') 

a 01 “ 

+ £ Ti Uk\i'k'\1'k')u ( jfcK I X "K " ) , 

a a 

(1.117) 


where we have defined 

T] (ZK |Z'k'|z"<") - 
a 1 1 


From the form of Eq. (1.117) we see that the angle change 
depends on the difference of displacements of atoms from 
the vertex atom. 

Direct computation using Eq. (1.115) gives 

cos 6 (o) - - 1/3 (1.119a) 


cot 0 (O) (XK \l'K' |Z"k')R (ZK |x'0 
|R(XK I I'K ' ) I 2 

CSC e Co) (ZK |X'K ' |X ,/ K")R a (XK |Z"k" ) 

|ft(ZK I X'K ") | |R(XK |Z # k") I 

» * 

(1.118) 


with 

9 ( 0 ) - 9^°^(X,0jZ + 6 1 ,1|X + 6 2 ,1) - 109.47122° . 

(1.119b) 

This angle is the tetrahedral angle and is the same for each 
of the twelve terms in Eq . (1.104). As an example, using 
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Eq. (1.117) we find 


A0(l,O|i+6 l|*+6 1) - 1_ £ 

v6 a 


( - 4 al- 5 O 2 +26 0[ 3 )u a (i - 0 l 1+6 l'« 

+ 4 4 al" t " 4 a2 +24 a3 * u a ^ » 0 U +i 2 * 1 J 

( 1 . 120 ) 


so that 


d 


dU a (l,0) 


2a r-> . 

3 f a3 


ar 


— (1 ' , 0 1 X '+ 6^, 1 1 1 '+6g, 1) 

X 


( - 6 ei- 6 e2 +2S s3 )u B (i ' 0|,+4 i’ 1) 


+ ( fi g2 + *32 4 ’^S3 ) U P ^ ^ I ^ +i, 2’ 


(1.121) 


This term only couples nearest neighbors. The first six 
terms in Eq . (1.7.04) are of this type. However the second 
six terms in Eq . (1.104) are different when differentiated 
as above. Consider 


A0U,i!x+6 5 ,o|x+x 6 ,) 


( - 6 al- 6 a2 +26 a3 >u a U ’ 1 l 6+ V 0) 


+ (6 ol +6 0 2 +24 a 3 )u a (i ’ 1 l i+1 6- 0) 


( 1 . 122 ) 


which gives upon differentiation 
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d 

3u a 


for! 

mxf. 


50 (i,l|i+6 5 ,0|i-f 6 ,1) 


-|^- s a l- t «2 +2i a3 ) 


(-6^-6^2+26^3)11^ ( i> ( 0 | i-63 , 1 ) 


+ (6 ai +6 02 +26 33 )u e U “ 6 5 +6 6 '°l i " 6 5 ‘ 1) 


- f ? (6 al +6 a 2 +26 a 3 ) 


^ 6 01 +6 0 2 +26 03 Ci-,0|X— 6^,1) 


+ ( - 6 01- & 02 +26 03 )u 0 U " 6 6 +6 5»°l i “ 6 6' 1) 


(1.123) 


This type of terra couples atoms that are second neighbors 
as well as atoms that are first neighbors , (i.e., the atom 
at ^(i-6g+6g, 0) is a second neighbor to the atom at 
ll(X,0) , and similarly for i?(i- A ,.+6g, 0) . ) Therefore we 
see explicitly how the "nearest neighbor" angle bending 
potential couples atoms that are first neighbors as well as 
atoms that are second neighbors. 

By successive evaluation of each term in Eq. (1.104) 
using Eqs . (1.117-1.118) and th^n differentiation of each 
term, applying Eq. (1.71) we can obtain the angle bending 
contribution to the dynamical matrix. This is a tedious 
exercise so we present the results here for reference pur- 


poses . 
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Dll } (0* 0 |it) - 


' 35< 008 S-7 2 + 5§ COS 5-T 3 1 

- A 008 5 ‘^4 " COS 5 ' ( W 

+ ll 008 < ^2"^4 ) "Bi 008 ( ^3 _ ^4 ) _ 

(1.124a) 


D^^O.OlJ) 


- 3TS 008 5 ' T 2 + A 008 £ ' ( W 

- W sln s ' <W + Ti r sln £ ' ( W 


aln 5,t 3 + W sin j,7 4 


(1.124b) 


D^3 > (0,0|it) 


- 3^5 cos S '*4 


cos £• (T 2 -T 3 ) 


2ic 

3M 


+ — sin £*T 0 - -?rrr sin k*T r 


+ ^sin K.(VV 'Tf sin * <VV 


(1.124c) 


D^ } (o,o|k) - d ^2 5 CO, 0 |lc) 


12 


(1 . 124d ) 


D^2 } (0» 0 i^) 


28c 

“ 3M cw ^ 

4c t* -*• . 4c 


F - & COS S ’ T 2 * #M cos 5 ' ”3 


+ 5§ COS k '^4 + 3M COS k ' ^2"^3^ 

- ^ cos k- <V ? 4 ) ' * C0S *’ ( W 


(1.124e) 
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D23 >C0 *0|^) 


- A 008 *3 + A cos ** ( W 

* TiT sia S' (W - ^ sin S. (£j-£ 4 ) 

+ S’VW** S’ ^ 4 


D 31 } ( 0, 0 |k) 

D^ } ( 0 , 0 |^) 

D33 } (0,0 (it) 


(1 . 124f ) 

D^ } (0, 0 |k) (1 . 124g ) 

(0, 0 (it) (1.124h) 

“3F + 3l COS k ’^2 A cos k ’^3 

’ A cos - A cos s ' ( V* 3 > 

/ A COS k ‘ ( W + A COS *’ ( VV * 

(1.1241) 


Eqs . (1.124d,g,h) were computed explicitly as a check on 
the algebra. They are also a consequence of the Hermit ian 
property of the dynamical matrix Eq. (1.42a), which must be 
obeyed for each contribution to the dynamical matrix. 

Similarly we have 


(0,1 |k) 


80 , . 12 ’ *2 lS-T 3 ^’*4, 

“ 3M^i + e + e + © ) 1 


Dl^ (0> 1 |S) 


(1.125a) 


A ik • ? 0 ik • ilc*?. 

40/4 ,2 3 4 X 

3Tm (1 + e " e “ e > » 


(1.125b) 
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♦ 


Dl3 } (©» 1 1^) 

D^^co^iJiT) 
D^CO.lJic) 
D^^O^lJk) 
D^® } (0,1 fit) 
D^2 } (0,l|k) 
D33 } ^0, 1 jk ) 


4o 

TKi 


(l - 


iic*T 
e * 


liT.^ 

e 


+ e 



(0,l|ic) 

Dll } (0»1 1^) 

1 Z't 2 ur.T 

3M 1 ” e +e " e 

Dl3 } (0» 1 ) 

D 23 } (0,1 | lc ) 

D ll } (0*1 1^) • 


(1.125c) 

(1.125d) 

(1 . 125e) 
(1.125f ) 
(1.125g) 
(1 .125h) 
(1.1251) 


Again, Eqs . (1.125d,g,h) were computed explicitly as a 

check on the algebra, they are also a consequence of Eq. 

(1.48), which must be satisfied for each contribution to 

the dynamical matrix. The t , in Eqs. (1.125-1.126) are 

give.: by Eq . (1.62a). Note also that our parameter a is 

23 

6Q in the notation of Herman. 



LONG-RANGE CONTRIBUTION TO THE DYNAMICAL MATRIX 


The possibility of comparatively long-range forces in 

diamond structure crystals has an interesting history. The 

first lattice dynamical investigation of the diamond struc- 

24 

turo was by Max Born. He considered the most general 
nearest neighbor interactions consistent with symmetry. 

The results for the elastic constants and Raman frequency 
in this model are 



C 


44 


(a-g 2 /q) 

a 



(2g-a) 

a 


8a 

M 


(1.126) 


Here a and 3 are the two nearest neighbor force con- 
stants. Since there are only two parameters in this model, 
this leads to a relation between the elastic constants 
known as the Born identity 


4C 11 (C 11 -C 44 ) 


(C 11 +C 12 ) 


and also the Raman frequency 



(1.127) 


(1.128) 
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I 
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When the first measurements of the elastic constants 

25 

for diamond become available. Born evaluated the left- 
hand side of Eq. (1.127) and obtained the value 1.10 in- 
stead of 1. This seemed rather good agreement, yet the 
ten per cent discrepancy caused Born to question whether 
second neighbor forces could improve matters. As it turned 

out, the value of 1.10 was fortuitous since later more 

26 

careful measurements showed the early value of C ^ to 

be more than a factor of three too large, the values of 

C, , and C AA also had large errors. Using the experi- 
ll 44 

07 

mental values for diamond, silicon and germanium, we 
evaluate Eqs. (1.127) and (1.128) and compare to experiment 


in Table 1. 
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TABLE 1 . 

Comparison of Corn identity Eq. (1.127) and Raman 
frequency Eq. (1.128) with experiment. Frequencies in 
radians/sec . 


Crystal 

*U< C U- C 44> 

(“■■)* 

Ui RA(expt . ) 

(C U^12 )Z 

Diamond 

1.49 

3 . 92xl0 14 

2 . 51xl0 14 

Silicon 

1.09 

1.24xl0 14 

9 . 79xl0 13 

Germanium 

0.995 

6.95xl0 13 

5.66xl0 13 
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As can be seen from Table 1, the agreement with experiment 
is the worst for diamond and the best for germanium. This 
shows the inadequacy of the nearest neighbor model to 
account for even the long wavelength vibrations of the 
crystal, though the agreement for silicon and germanium is 

respectable considering the simplicity of the model. 

28 

Helen Smith, collaborating with Born, worked out the 
theory for the case of first and second neighbor interac- 
tions. Her second neighbor force constant matrices are not 

29 

the most general allowed by symmetry and neglect the 
antisymmetric second neighbor force constant. Using a 
three parameter version of the model developed by Smith, 

o 0 

Hsieh fit the three elastic constants of silicon and 
germanium and computed the specific heat of these crystals. 
Hsieh found that the calculated specific heat values were 
well below the experimental values, even at temperatures 
as low as 60K , the discrepancy increasing with increas- 
ing temperature. This was an indication that the normal 
mode frequencies in the rest of the Brillouin zone were 
considerably smaller than those expected from the first 
and second neighbor model. 

When phonon dispersion curves became available from 
neutron scattering studies, detailed comparison between 
lattice model predictions and experiment was possible for 
short wavelength modes. Herman presented an analysis of 
diamond structure crystals using general forces between 
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first through sixth neighbors. He shoved that it was 
necessary to include interactions to at least fifth neigh- 
bors in order to fit the three elastic constants, the Raman 
frequency and the LA and TA frequencies at the [HI] 
zone boundary of germanium. 

In attempts to reduce the number of parameters re- 
quired to fit the phonon dispersion curve data, several 

32 

lattice dynamical models have been Introduced. Lax in- 
cluded displacement induced quadrupole-quadrupole inter- 
actions between atoms. This was an attempt to include long 

range interactions within the adiabatic approximation. The 
33 

shell model also includes long-range electrostatic in- 
teractions, but introduces electron coordinates which are 

34 

then eliminated adiabatically . The bond-charge model is 
similar to the shell model but employs point charges be- 
tween the ions that are allowed to move adiabatically. 

35 

Finally there are the valence force models which describe 
the forces between atoms in terms of bond stretching, bond 
angle bending and combinations of the two. 

In each of the above models it was found necessary to 
include interactions of some type to at least fifth neigh- 
bors in order to fit the phonon dispersion curves and the 
elastic constants. Each of the models has attendant 
strengths and weaknesses. Since the author wanted a har- 
monic model to use for computing anharmonic properties, 
none of the above models were considered satisfactory for 
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this purpose. Therefore we have investigated a model sug- 

36 

gested by the work cf Lax, but with very different con- 
sequences than his quadrupole model. 



NONLOCAL DIPOLE MODEL 


The presence of a lattice displacement can distort the 
electron charge density in the vicinity oi a displaced ion. 
This induced charge density can then Interact electro- 
statically with a charge distortion produced by another 
displaced atom, giving rise to long-range interactions. 

How important this effect is should intuitively depend on 
how big a charge distortion is induced by a given lattice 
displacement, and how effectively this charge distortion 
is screened. 

To put this concept on a firmer foundation, in Figure 

6 is plotted the valence charge density of diamond and 

silicon, as computed ab initio from the electronic theory 
3 7 

of solids. This is for the lattice with all atoms in the 
rest positions. In Figure 7 is the induced charge density 
for silicon due to the presence of a [100] transverse zone 
boundary phonon. The induced charge density is the dif- 
ference in charge densities between the distorted and un- 
distorted lattices. Note that there is a striking dipole 
contribution to the induced charge density, as well as con- 
tributions from higher multipoles. Since we only want to 
model the longest range contribution to the forces, we will 
examine the dipole contribution to the induced charge den- 
sity in detail. 

We begin by writing the dipole moment induced about 
atom (4K) due to lattice displacements as 
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Figure 6. 
silicon. Units 


Valence charge density for diamond 
are (e/^ 0 ) • (Courtesy M.T. Yin, 


and 

Ref. 37) 
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Valence charge density { 1 10 plane) 



Figure 7. [ 110 ] plane charge density Induced in sill* 

con by the presence of a TAX phonon. Units are (e/fl Q ) . 

0 

Amplitude of the phonon is .077A . Atoms not matched cor- 
responds to taking the difference in charge densities with 
and without the phonon using the same origin of coordinates. 
Atoms matched corresponds to taking the difference in 
charge densities with the origin of coordinates for the 
phonon present charge density shifted so that the atoms 
shown are in the original undistorted positions. (Courtesy 


M.T. Yin, Ref. 37) 
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(a) TA(X) 

(Atoms not matched) 



(b) TAM 

(Atoms matched) 
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> 


I 


p (iK> - E p -(4K|l'K # )u fl U'K # ) 
a X'tt'p uP p 


(1.129) 


38 

This form has been mentioned by Lax to justify his 

quadrupolo model. Although the form of Eq. (1.129) Is 

39 

similar to the work of Minnick, our interaction energy is 
not the same as his and leads to rather different conse- 
quences. Since Eq. (1.129) had not been applied to the 
diamond structure in the manner that we apply it, it was 
not known what effect it ould have on the phonon spectrum. 
Therefore we have not considered polarizability effects, 
which in our approach requires considerable additional 
effort to include. 

The total dipole moment of the crystal may be written 
as 


2 p ( X K ) , 

u a 


(1.130) 


so that. 


p - E E p fl (XK |i'K')u, (X'O . (1.131) 

a U i-K'S & ? 


Defining 


M ae< rK') - z p ae UK|ro , 


(1.132) 


we have 


E 

X'k '0 


(X'K')Ug U'k:') . 


(1.133) 
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The transformation property of the total electric moment 
under the operation. (Sjv(s) +R(m)} (see Eq. (1.34)) leads 
to 40 

t V LK) * S uo S uS M aP <1K) ' (1 ' 134) 

OLP 

The special case of a pure lattice translation by it(m) 
gives 

M (l+m,K) - M^U,K) - M^,(0,K) , (1.135) 

so that M , (IK) does not depend on the unit cell index. 

Hi' 

Since an arbitrary uniform translation of the crystal as a 
whole must not change the total dipole moment, we obtain 
the requirement that 

2 M a (i ' o - 0 . (1.136) 

t'K' 

Using Eqs. (1.135) and (1.136) we also have 

£ M a (0, K) = 0 . (1.137) 

K 

This is the statement that the primitive unit cell is 
neutral if M a ^(0,K) is thought of as some type of effec- 
tive charge. Applying the inversion operation Eq. (1.36) 
yields 


M ag U,°) - M ap (-i,l) . 


(1.138) 
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With the addition of lattice translation invariance^ Eq. 
(1.135), we obtain 

M og (0,0) - M a0 (O,l) , (1.139) 

so that Eq. (1.137) becomes 

M a0 (O,O) - 0 , (1.140a) 

or 

M ag U,K) - 0 . (1.140b) 

Eq. (1.140b) expresses the fact that no linear term in the 

total electric moment exsists in diamond structure crystals. 

This has the consequence that the infrared absorption 

spectrum must be explained by two and higher phonon pro- 
41 

cesses, although the presence of impurities and surfaces 
can break the translation invariance and allow a first 
order moment to exsist. 

From Eqs . (1.140) and (1.132) we obtain for diamond 
structure the requirement that 

Z - 0 . (1.141) 

u 

Note that Eq. (1.141) does not require that p^(XK) be 

zero in the bulk crystal. This v?ould only result if we 

had used a strictly local form for p^ 0 (Xk | £ ' k ' ) , as used 

42 

in the work of Trullinger. 
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The transformation property of the dipole moment 

-► ^ 43 

p^CiK) under the operation [S | v(S) + R(m) } leads to 


p^(LK|L'K') - S . (1.142) 

The special case of a pure lattice translation by it(m) 
gives 


P a(3 (X+m,K |i'+m,K') - p ag (XK|i'K') , (1.143) 


so that in particular 

p ae l X#K “ p a0 (0 ' K ■ P a g U-l't* |o,k ') . 

(1.144) 


Thus we see that the dement coefficients transform in che 
same manner as the second order force constants Eq. (1.35). 
However, the moment coefficients (IK j i *k ' ) need not 

obey the "flip” symmetry property Eq. (1.26) that is obeyed 
by the second order force constants. 

Consider now the moment coefficients appropriate to 
the case of nearest neighbor nonlocality. This is the 
simplest type of nonlocality and is the model we will in- 
vestigate. Applyinr 'q. (1.129) to this case we have 


P a U,0) 


£ p ttg U,0|£,0)u a U,0) 

0 

4 

+ £ £ p „U,0|4+6 1 )u r U+6.,1) 

0 i~l a P 1 P 1 


(1.145a) 



8 

+ E E P„ fl U,l|X+&,,0)u fl U+6.,0) , (1.145b) 

8 i-4 “ p p 1 

where the 6^ are given by Eqs . (1.62). Translation in- 
variance plus inversion symmetry analogous to Eqs. (1.40) 
yields 


p eo (/,o!i,o) - p oB <o,o|o,o) - p aB a,i|z,i> u.i46a) 


and 


p p U,l|i+6 i ,0) - p (0,01-6^1) 


(1.146b) 


Using Eqs. (1.146) and (1.62b) we obtain 


P a (A,0) 


E P ag (0,0|0,0)u g (1,0) 
0 


+ E E p fi (0,0|6,,l)u R (X+6.,1) (1.147a) 

8 i*l ap i P i 


P (X,l) - 2 P ag (0,0|0,0)u g u,l) 
0 


+ E E p p (0,0|6 i ,l)u g U-6 i ,0) . (1.147b) 
8 i~l 


Applying the requirement Eq. (1.141) and inversion sym- 
metry analogous to Eq. (1.40c) yields 


p aj3 (0,0|0,0) 


- E p (0,0|6.,1) , (1.148) 

i-1 ap 1 


so that 
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p (1,0) - £ £ P a g{0,0|6 i ,l)(Ug(i+6 i ,l)-Ug(i,0)) 

3 i"l 


(1.149a) 


p (1,1) - £ E P ag (0,0|6 1 ,l)(u p (i-6 i ,0)-u p U,l)) . 
8 i-1 


(1.149b) 


Using Eqs. (1.142) and (1.68), we can immediately write 
down the four matrices 


p(0, 0 1 (* 1 ,1) - /p 1 


p(0,0|6-,l) - 


P 1 P 2 ” P 2) 


p 2 P 1 P 2 


V p 2 p 2 P l/ 


P 2 P 1 ” P 2 
\" p 2 “ P 2 P V 


P(0,0|6 3 ,l) 


Pl -p 2 ” p 2' 
” P 2 P 1 P 2 
V- p 2 P 2 P l- 


P ( 0, G j i 1 ) ■ / P 1 ” P 2 P 2 

“ p 2 P 1 ” P 2 

{P 2 ~ P 2 P l' 
(1.150) 


Here p^ and Pg are parameters to be determined and 

have the dimensions of charge. The matrices Eq. (1.150) 

44 

agree with those of Lax. Using the explicit form of the 
matrices we obtain 

£ p ae <°<°l 6 i> 1) - 4 Vce - (1 - 151 *> 

so that Eq. (1,148) becomes 

p a .(0,0|0,0) - - 4 Pl 6 ae . 


(1.151b) 



At this point we note that Eqs. (1.149) satisfy in- 
finitesimal translation invariance explicitly, though this 

condition was not imposed on p (IK) in our derivation. 

o 

Only the total moment Eq. (1.130) was required to be trans- 
lation invariant since we are dealing with a neutral cry- 
stal. Lattice translation invariance plus inversion sym- 
metry then gave Eq. (1.141), which is not the same condi- 
tion obtained by imposing infinitesimal translation 
invariance on the moments P a (lK) . In fact, the defini- 
tion of p (IK) requires a volume to be associated with 
a 

each atom (IK) . If this volume is neutral, then the 
dipole moment is independent of the origin of coordinates 
and is translation invariant. If the volume is not neutral, 
then the dipole moment does depend on the origin of coordi- 
nates and is not translation invariant. The fact that 
Eqs. (1.49) are translation invariant without imposing this 
condition is a reflection of the fact that the first neigh- 
bor moment coefficients are symmetric and this together 
with the inversion symmetry shows that 

P aS (0,0 f ^i* 13 " Ppa^i ' 1 f°,°) • (1.152) 

Thus we see that for first neighbors, the "flip” symmetry 
condition Eq. (1.26) is obeyed by the moment coefficients. 
For second neighbors this condition need not be true. When 

v 

the "flip" symmetry property is obeyed, Eq. (1.141) is then 
the same condition obtained by imposing infinitesimal 
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translation Invariance on the moments p (XK) . 

at 

One can explicitly verify that the total dipole moment 
due to Eqs. (1.149) is always zero, as we have seen it must 
be from Eqs. (1.133) and (1.140b). Although the total 
dipole moment is zero the array of dipoles can have a total 
quadrupole moment as shown explicitly in Appendix C. 

Having determined the form of the moment coefficients, 
we now proceed to calculate the dipole contribution to the 
second order force constants. In Appendix B it is shown 
that for silicon we can write the interaction energy for a 
pair of dipoles as 


wuK|x'o-i £ n UK|rOp uop-u'o 

€ a,g ap a p 


(1.153) 


with 


n a 0 UK|ro 


ap 


R(XK Z'K') 


3R (XK | X'k ')R a (XK X' k') 

- H 

1 3 


|ft(XK |X 'k ') 


(1.154) 


Here e. is the static dielectric constant. In writing 
Eq. (1.153) we have neglected short-range terms discussed 
in Appendix B. These terms are compensated for by the 
short-range contributions to the dynamical matrix. 

The total dipole-dipole interaction energy is 

4 dd - \ 2 £' W(XK | X'k') 

2 XK X'k ' 


9 


(1.155) 
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where the prime on the sum means exclude terms with 
(iK) - ( I'k ') . Using Eq. (1.17) we have 


^l K i I A 2 K 2^ " \ ^ ^ ' 

<*P 1 1 2 2 2 i Kn x * 


3 3 '4 4 


J 2 


3 W <^3 K 3'V4 ) 


dU a (ii“i)du 0 (i 2 K 2 ) 


Jo 

(1.156) 


where the derivatives are evaluated with all atoms in the 
rest positions. Using Eqs. (1.149, 1.153 and 1.154) it is 
a straightforward excercise to evaluate the derivatives in 
Eq. (1.156). We obtain the results 

dd 16p l 

*ae<V°!V°> “ -T in ae< 4 i»°l A 2' 0)<1 " 6 i 1 ,i 2 ) 

, 4 

+7 2 2 p (0,0|6 l)p (0,0|6 ,1) X 

€ i,j-l p J 


X ^ l / ( V 6 i^ 1 l i 2 +6 j' 1)(1 " 6 

4 


WV 6 .1 


-7 2 Lp (0,0|5 l)p (0,0|6 1) 

€ i,j-l p.,1/ p J 


X 


x [^<i 1 .0| V 6 J> 1)+ ‘V U 2' 0|i l +6 i' l) ] ' (1 ' 157a> 


and 
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_ 16p7 

•a? ( ' t l- 0 IV 1) * -T in a e (i l'°l 1 2 ' 1) 

+ « i j-1 ^a/ 0 >°N> 1) v <0 ' ol v 1,fi 5l' < v s l' 1 !v 6 j> 0) 

J p <8i »'. 0 l‘t* 1 > n |1 ii<V 6 i« 0 lV 0>(1 -V‘i* 1 2 > 

‘ "T £ S V <0 ' 0|s i’ 1)n li a U l +4 l' 0 l 1 2' 0>a -V 6 i' 1 2 ) ' 


The Kronecker delta symbols in Eqs. (1.157) mean omit terms 
when the Kronecker delta is one. Defining 

d£?(kk'|S> - (1.158) 

X 


we find 


»£«>.o|E) * n ae( 1 - 0 | 0 ’ 0 ,e ' ik "' U) (i' 6 i.o) 

+ i S Sf *»*„«)« ; (i t 0 |°. 0 ).- 1 * ,f( 1 , (i-» lf0 ) 

* \i,V l 

M> * 

-^LEf (it)a B u,o|o,i)e +ik, ^ u) , a 

€M , aw- M-P 
M- * 


and 
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D a8 (0 ' 1, ^ ) " 'TST 2 

i 

+ 7B ^ S F aii (S)F ei/ (f ) n iil/ (x,i|o,o)e- lI,S(<) 

- F (m (S)«\ rf «.0|0.0).- lK, *< i >(1.6 lj0 ) 

- TirSE F 3 /^%(^O|0,0)e- ie,S(i >(l-6 i(O ) 


where we have defined 


F “ s<e) ' J 1 p a s ( 0 ’ 0 , 6 i’ 1,eii ‘' 7t • 

Decomposing F ^ (k) as 

F aS < S) * [PiS(IU ae +P 2 S a| 3 (k) 1 


we have 


ik*T 9 lit* 5_ ik» t 
S(k) - (1+e J +e J +e 4 ) 


and 


Sab 






9 


(1.159b) 


(1.160) 


(1.161) 


(1.162) 


(1.163) 


with 
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’12 


*13 


’23 


iii'T, iJc • t- ik*> 6 a 
(1+e 2 -*e 3 -e 4 ) 

i£* 6« lk' i£»T 4 
(i-e 2 -e 3 +e 4 ) 

iis* T Q lie*?- ik’ t, 
(1-e 2 +e -e 4 ) . 


(1.164a) 

(1.164b) 

(1.164c) 


The dipole contribution to the dynamical matrix can be 
written in a more compact form as follows 


D^(0,0|k) 


4irp l _11^. 4 * P 1 P 2 *12,^ 

rSKT T <*0 (k) + eMft T oP (k) 


47r Po 22 -+ 
+ THTT" T aS <k) 


(1.165a) 


D o6 <°’l !^) 


. llfi V 11 ^) + 4TP )- P 2 v 1 ? (k) 

TmJT v « 0 U5; eMT. a0 v ' 


4irp 


2 y2? 


€MC o a0 


(k) . 


(1.165b) 


The dimensionless matrices T and V are defined as 


T “(S) 


^~l [l6+S*(k)S(k) ] Z' (X,0|0,0)e 


[- 4s*(it) (i,o|o,i)e -ik, ^ (X) 

JL 

- 4s(iT) 2n a3 u,o|o,i)e +lk * 5(x) 


-iiMtu) 


(1.166a) 


82 


n 


- £ 


stfjEs^dfjE'o i i< „<i,o|o,o).-‘ s, *<*> 


l M.Of 


s*(?) Es^ (?) E' u,o|o,o)e” ik,R(x) 

M* i 

“ 4 Ss ^^) 2 n ua u,c|o,i)e“ lk,R(x) 


" 4 SS au (i ‘ ) 2 :n u 0 U ^ o l O ' 1)e+iJC ' RU) 

u * l J 


(1.166b) 


T J <k) 


n 


’ 47 ^ S a 4 (k)S g^ ( ^ ) ^‘V Ul0f0 » 0)e ’ ik,RU) 


(1.166c) 


V 1 ? 

ap 


(k) - 


n 

17 ] 


16 Ef^ U,OjO,l)e" ik ’ R(X) 


+ s 2 (?)£n ag u,i|0,0)e“ i ^ ,R(x) 

JL 

r BstoS'fi^ U,0|0,0)e‘ tS,irU) 


(1.167a) 


Q 


V 12 (?) . Ij> 
4tt 


S(if)2s p>i ()f)£n ii(;[ (i ( l| 0 ,0)e- ie ' S ' u) 1 


+ S(k) 2 s (if)£n (i, l | 0 , 0 )e’ ik,RU) 

M. Z 


- ■‘ 2s mi < k >2'n e (i.o|o,o)e* - 

M. I 


itu) 


L- 42s j 3 4 (^£' n ^(i)°| 0,0)e“ i ^ , ^ (X) 


(1.167b) 
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V^k) - ^s w <ic)s ei/ (C) 2a u „(x,i|o o)e- i^ • S(^, . 

(1.167c) 

In Eqs . (1.166) and (1.167) there are three long- 
range buds to perform, two of which are related, namely 




«<?’<*> 




^2 £' n ag u,0|0.0)e' 1 ' t '^ <<) 

X 

7 ? 2 Q mJ <Z,l|O,0)e' lk,RU> 
t 


(1.168a) 

(1.168b) 

(1.168c) 


These sums are slowly convergent for numerical work. In 
Appendix D we discuss the Ewald method of transforming 
these sums to rapidly convergent form. The Ewald method 
also has the advantage of allowing explicit separation of 
the nonanalytic behaviour of Eqs. (1.168) which enables one 
to show that the dipole contribution to the dynamical mat- 
rix is analytic at lc - 0 „ even though the separate con- 
tributions to it are not. 


r 


PHONON FREQUENCIES ALONG SYMMETRY DIRECTIONS 

In this section we present expressions for the dynami- 
cal matrix and for phonon frequencies along the [100], 

[110] and [111] directions. We also give expressions for 
the elastic constants using the model developed in the 
previous sections. 

[100] Direction 

Along the [100] direction the wavevector lc is 

k - ke 1 . (1.169) 


Direct computation using the model developed in previous 
sections allows one to write the dynamical matrix along 
[100] as 



(1.170) 


where and Ag are real and we have defined 

A x - A + D^(0,0|k) (1.171a) 

A 2 - B + D 22 (0» 0 |k) (1.171b) 
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z_ 
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A 3 ■ C + Djj(0,l|k) 


a 4 - D + D^(0,l|k) 


A 5 - E + ^23 ( 0 » 1 1 k ) 


(1.171c) 
(1 . 171d) 
(1.171e) 


and 

A - 


C = 

D - 


4a4-16u sin 2 ^ + 8p. *+4\ * 4-4X *sin 2 4- 12c - 4p cos ^p 


(1.172a) 


B - £ 


4ar4-8 ( m- 4-X )s in 2 ^p 4- (8u +4\ )+4p, sin^^p 


2 ka 


1 

M 

1 

M 


26ct 0 ka 

+ —j" + 2c cos ~y 


/n . ' , 16a w1 , -ika/2 x ika/2, -Ikav 

( 2q'4-4|j. +-^—)(l4-e )4-2X (e +e ) 


(1.172b) 

(1.172c) 


2or4-2X. * 4-i|^) (l+e" ika/2 )+2^'(l4-e ika/2 + 


-ika/2 , „ -ika x 
e 4-e ) 


(1 . 172d) 


E - t. 


(.26+46 ' + ^ (1 . e - lka/ 2 )+2t , ' (e ika/2 -e- lka ) 


(1 . 172e) 


We now proceed to find the eigenvalues of Eq. (1.170) 

45 

using the eigenvectors of Lax. It is important to note 
that his eigenvectors correspond to the second form of the 
dynamical matrix ^’lich uses displacements of the form 


x _ V"> ilc • ft( XK )-iurt , 

' r— “ © 


u * <iK ' Vm 


(1.173) 
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and differs by a phase factor from Eq. (1.19). Thus we 
have 


w (K) - V (K)e ik '^ (K) , 


(1.174) 


where W (k) is defined by Eq. (1.19). The V (k) are 
Ot ” 

given by Lax for the [100], [110] and [111] directions. 

Consider first the longitudinal modes. Noting that 
along [100] 


( ik*K(l) . ika/4 


(1.175) 


we have for the longitudinal acoustic (LA) mode 


W (0) - b6 , 
o al 


w \ ika/4. , 

V 1 ’ - e b6 ol . 


(1.176a) 

(1.176b) 


Substituting Eqs . (1.176) into Eq . (1.20) and using Eq . 
(1.170) we obtain 


2, ika/4.. 

u> b - (A^^ +A 3 e )b 


2, , . t .* -ika/4.. 

a: b - (A 1 fA^e )b 


(1.177a) 

(1.177b) 


For this to be a solution we must have that 


A 3 e 


ika/4 


purely real. 


(1.178) 
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This Is in fact the case so that 


<*[1001 - [A 1+ Re(e lW4 A 3 )l 


(1.179) 


After some algebra we obtain 


2 (8a+16u'+64c/3) A _2 ka 

“LA(IOO) ‘ ■ i “ a M al " X 

+ 06^/32 sl „2 + 8^1 aln 2 3|a 

2 


+ sin" ^ + ^i[TjJ(k)+Re(s ika/4 vJJ(k))] , 
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M 


2 ka 


(1.180) 


and we have used the facts that 

T ll (k) " T ll (k) " V ll (k) " V ll (k) " 0 U-181) 

along [100] and 


Tn (It) - purely real. 


(1.182) 


We can obtain an expression for the elastic constant Cn 
from Eq. (1.180) by using the fact that for ka « 1 


^[ 100 ] 



higher order 
in ka 


(1.183) 


aCn - (a + 4c + 8ii + 2u ' + 91 ' + 8\" ) . 


(1.184) 


so that 
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Note that the dipole terms are absent from Eq. (1.184). 

This is an important point and is a consequence of the non- 
local nature of the dipoles. In order to see this, we only 
need show that along [100] 




+ Re(e 


ika/4 v ll 


V“UO) 


- 0 


(1.185a) 


k-0 


and 

- 0 . (1.185b) 

k-0 

This is in fact the case as shown in Appendix E. 

Consider next the longitudinal optic (L0) a' mode, 
for this mode 

W (0) - b6 , (1.186a) 

a al 

W (1) - - e lka/4 b6 . . (1.186b) 

a al 

Substituting Eqs . (1.186) into Eq . (1.20) and using Eq. 
(1.170) we obtain 


a 

ai^ 


T^(k) +Re(e ika/4 V^(k)) 


2 W , . . ika/4v. 

oj b - (A^ - A^e )b 


(1.187a) 


2, / a .* -ika/4 N . 

cu b - (A^-A^e )b 


(1.187b) 


using Eq. (1.178) we have 
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^0(100] • [A l -Re(e U “ /4 A3> 1 ' (1 ' 188) 

After some algebra we obtain 

2 _ (4a+8u '+32o/3) „ ka x , 16m, _^_2 ka 

“LOtlOOi if i(l+co«- r )+-jf»m -J- 

4\' .. , 3ka N , 4\ " , 2 ka , 4 a n ___ ka v 
+ _ (l + c °s _) +-JT sin -y + 3M a -cos^) 

2 

+ 7sn itT ii (k) “ Re(eika/4v ii (]k)) ] • t 1 - 180 ) 

o 

I becomes the Raman frequency , denoted 

, so that 

- (8a + 64a/3 + 16u ' + 8\ ' ) /M . (1.190) 

In obtaining Eq. (1.190) we have made use of the fact that 

T^(k-O) - Re V^(k-O) - 0 . (1.191) 

Note that the dipoles do not contribute to the Raman fre- 
quency, which is again a consequence of the nonlocal nature 

of the dipoles. This is in marked contrast to the local 

46 

quadrupoles used by Lax which do contribute to both C 11 
and . 

Furthermore, at the X point of the Brillouin zone 
boundary k * , so that we find from Eqs. (1.180) and 

a 

(1.189) 


At 


K - U 


^ 0(100 
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u£ ox - “ |<4a^0c/3+16u+«u'+4X') 

4 irp? 1 1 

+ «■£ *£<*> <1 ' l92 * ) 

where 

T^;(X) - 5.51804915 . (1.192b) 

ii 

In obtaining Eq. (1.192a) we have used the fact that 
V**(X) - 0 . 

Now we proceed to the transverse modes. These modes 
are labeled and are doubly degenerate. One eigen- 

vector is 


V 0) ■ a(6 0 t 2 - 6 o 3 > (1.193a) 

*„<!> ' e lka/4 b(6 a2 " 6 a 3 ) • (1.193b) 


Substituting Eqs . (1.193) into Eq. (1.20) and using Eq. 
(1.170) we obtain the secular equation for the transverse 
modes 


^[lOOji " A 2 * t (A 4“ A 5 )e 


ika/4, A *_ A *^-ka/4j2 


(A 4 - A 5 )e 


(1.194) 


We can rewrite Eq. (1.194) as 
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<4[100]* • B + D M (0 '°l k) * 

[ Re (De lk,/4 >+Re <e Uk/4 D^ ( 0 , 1 1 k ) ) ] 2 1 * 
+[lm(Ee lka/4 ) + In(e ika/4 D^(0,l|k))] 2 ) (1.195) 

where we note that 

De lka./4 _ _ i |^(4a+4pi '+4X '+32a/3 )cos ^ + 4u ' cos “pj 

(1.196a) 

£e ika/4 _ i^(i6o/3^«6'-4@)sin ~ + 4i/'sin^ 

(1.196b) 


and 

D^^O^IkJe^ 74 - Re(e ika/4 D^(0,l|k)) (1.197a) 

D 23 (O’ 1 |k)e ika/4 - lIm(e ika/4 D^(0,l|k)) . (1.197b) 


To obtain an expression for the elastic constant C 44 
we use the fact that for ka « 1 


2 

^AflOO] 



higher order 
in ka 


) 


(1.198) 


so that 


I 
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aC^ • o + + 4u + 4X + X' + 10u' + 8>i # 


(4o/3-0 +26 '+3p' ) 2 
” (a+da/5+iii'+x p ) " * 


(1.199) 


In obtaining Eq. (1.199) we have used the tact that the 
lower sign in Eq. (1.195) corresponds to the acoustic 
branch. Note again that the dipole terms do not contri- 
bute to C 4 4 . 

At the X point of the Brillouin zone boundary we 
find from Eq. (1.195) that 

U^AX • ^(4a-4g+8|i+8X+8u '+4X ' ~4v '+8 6 '-*-12 a) 

4v(p--p«) ii 

+ wrr — ( t 22 ( x ) + v 23 < x)) (1 - 200) 

o 

and 


2 

^TOX 


( 4of+40 -+8p, +8 X +8n '+4X '+4v '-86 '+4c/3) 

M 


4tt 


(P 1* P 2> ,.11 


.11 


€ Mu 


( T 22 (X) “ V 23 (X)) 


( 1 . 201 ) 


where 


T 22 (X) " "2 . 75902458 (1.202a) 

v£*(X) - -13.5224054 . (1.202b) 


In obtaining Eqs. (1.200) and (1.201) we have used the 


relations 


(1.203a) 


vJJ(x) - vi^(x) 

T 22 (X) - T^W) (1.203b) 

T 22 (X) " “ 2V 23 (X) (1.203c) 

vJ 2 (X) - “^T ^2 (X) . (1.203d) 

The dipoles strongly affect the zone boundary modes. They 
push the transverse optic (TO) mode up and lower the 
frequency of the transverse acoustic (TA) mode. 


[110] Direction 

Along the f 110 ] direction the wavevector k is 


k 


k „ 

S2 ie 1 


+ e 2 ) 


(1.204) 


where k - |k| . Direct computation using the model 
developed in previous sections allows one to write the 
dynamical matrix along [110] as 



(1.205) 
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r 


94 


where , Bg and B 4 are real, B 3 la purely 
and we have defined 


B 

B 

B 

B 

B 

B 

B 

B 


1 

a. 

3 

4 

5 

6 

7 

8 


A + D 
B + D 
C + D 
D + D 
E + D 
F + D 
G + D 
H + D 


dd 

11 

dd 

12 

dd 

13 

dd 

33 

dd 

11 

dd 

12 

dd 

13 

dd 

33 


o 

% 

o 

'w' 

|k) 

(0,0 

|k) 

(0,0 

|k) 

(0,0 

|k) 

(0,1 

jk) 

(0,1 

| k ) 

(0,1 

(k) 

(0,1 

|k) 


and 


A - - 


4o£+ 8 )sin 


2 ka 

73 


+ 4j* sin 


2 _ka_ 

2v^ 


8u' 


, , „ » . 2 ka „ _ ka 

+ 4(n +X )sin - 9a r 2a cos 


- a/3 cos yj 


B “ ^j(4y + 2a/3)sin 


2 ka 

172 


i 4a ka 2a ka v 

' M ( T sin i7f-T sin 75 ) 


imaginary 

(1.206a) 
(1.206b) 
(1.206c) 
(1 . 206d ) 
(1 . 206e ) 
(1.206f ) 

(1 . 206g ) 
(1 . 206h ) 


J (1 .207a) 
(1.207b) 


C 


(1.207c) 
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if 2 ka 2 k* «i 

D - i 4a + 16 ^ Bin* 6 ^ + 4X *in ^75 

, , , 2 fc* 32o 

+ 8*' + 4\ + 8 »' .in 45 + — 

4a 4a J£*L 

L+ T 008 yz ’ T C ° 8 2/Z J 

E . - Ir a ^a/3)a + 2e- lkl/2V5 «>- lka/ ' /! ) 


(1 , 2u7d ) 


+ u ' (2+ 2e- Ut » /2 ^ 5 + 2 e - lka/ ' /5 + e lkl/2 ^«- 31kl ‘ /i;5 


, ,1 ^-ika/VZ ika/2yZ -3ika/2«/Z, 


+ X (1+e 


F - (-e+4o/3)(l-f«“ ika/>/? -2e“ ika/2 ^) 

+ 26 4 (l + e- ika/ ^e ika/2 ^-o- 3ika/2 ^) 


(1 . 207e ) 


- v * (l-2e” ika/2 V^ +€ “ ika/ '^) 


(1 . 207 f ) 


i[(-0+4a/3)(l-e" ika/ ^)+6'(e ika/2 ^-e” 3ika/2,/J ) 


, ika/2^5 -ikaA/Z -3ixa/2,/Z. 
1 / (1+e -e -e J 


~ _ i f(a+X '+8o/3 ) (l--e“ ika/ ^+2e'’ lka ' /2v ^) 


(1 . 207g ) 


+ 2n'(l+e“ ika/ ' / ^+« ika/2 ^ 5 +e" 3ika/2 ^)J . (1 . 207h ) 


Now consider the modes whose eigenvectors are deter- 
mined by symmetry alone. Noting that along [110] 
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^ika/zyaf 

w " 6 


(1.208) 


we have for the £. mode which la TA _ 

4 XZ 


V°> - >(6 al -6 a2 ) 


"a a) * « lka/2v,2b < 6 o 1 - t a2 ) • 


(1.209a) 

(l.C09b) 


Here TA denotes that the £. mode la transverse 
xz 4 

acoustic and polarized perpendicular to the z axis as 
well as to the diiection of propagation. Substituting Eqs. 
(1.209) into Eq. (1.20) and using Eq. (1.205) we obtain 


A - [B 1 -B 2 +e ika/2 ^(B 5 -B 6 ) ]b 


For this to be a solution we must have that 


(1.210a) 


A - [B 1 -B 2 +e“ lka/,A/J (B5-Bg) ]b . (1.210b) 


e ika/2yZ (B ^_ B ^ } _ pure ly real 


(1.211) 


This is in fact the case so that 


uj£ - Dj-Bg + Re[e ika/ 2 v/ 2 (B 5 -B 6 ) ] 

4 


( 1 . 212 ) 


After some algebra we obtain 
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" |jj(4of 48-t8n+8X+8n 4 p *+4\ '+&b *+12o)sin 4^5 

+ |j(4n-4t'+4^ *+4X *-8 6 *+4\k # +4X " )ain 2 yjj? 

+ DjJ(0,0|k)-Dj2(0,0|k)+Re(e lka/2v/5 DjJ(0,l|k)) 

- Re(e lka/2v ^D^(0,l|k)) . (1.213) 


Note that D°°(0,0|k) and 0 1 k ) are purely real 

along [110], We can obtain an expression for the elastic 

constant C ^ by using the fact that the £ 4 mode is the 
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"slow" TA mode in the elastic region so that for 
ka « 1 


2 

“s. 


(Cn-Ci 2 > 2 I higher orders 


Sp 


\ in ka 


Using Eq. (1.213) and the expression for , 

wo find that 


. . (1.214) 

Eq. (1.184), 


aC 12 " (20-a-2c-i-8i'“4^-4\-lOM, '+2v '-X '+126 '-8p, " ) . 

(1.215) 


From the expression for the bulk modulus B 


B 


3 (C 11 +2C 12 ) 


(1.216) 


we obtain using Eqs . (1.184) and (1.215) that 



» 


98 


B - ^[40-a+16l/-fc'X-18vi%4v'+7X # +24t'-^8\'-164'j . 

(1.217) 

Note that the dipole terms do not contribute to C^ 2 and 
that the bulk modulus is independent of the angle bending 
interaction. We can see in a simple way why the dipoles 
do not affect C 12 • Since the dipoles do not affect , 

if they do not affect the bulk modulus, then by Eq. (1.216) 
we see that they will not affect C^ 2 . The reason that 
the dipoles do not contribute to the bulk modulus follows 
from the fact that a homogeneous isotropic deformation in- 
duces no dipole moments and thus no dipole-dipole inter- 
action energy for this type of deformation. 

Next consider the mode which is TO . For this 

2 iz 

mode 


V°> 


W (11 

a 


b(t al- s «2> 


ika 2 -'\(-6 +6 ) 

oi a2 


(1.218a) 

(1.218b) 


Substituting Eqs . (1.218) into Eq . (1.20) and using Eq . 
(1 . 205 ) we obtain 


u- ,2 b - [B^B^e 1 **' ^(Bg-B^ ]b 


(1.219a) 


u) 2 b - i B i" B 2 + e ~ ika ' 2 ^(Bg-B* ) ]b . (1.219b) 


Using Eq . (1.211) we have 
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u 2 - Bj-Bj + He[e lk * /2v/5 (B 6 -B 5 ) ) . (1.220) 

Z 2 

After some algebra we obtain 

"" — 

(6a+10n *+2v '+4X '+34a/3 ) 

+ (2a+4u'+46'+2X'+10o)cos ^ 

lc& 

•+• (2ii +2X — 2v -46 )cos 

2 ka , . 2 ka 

+ ( 40 +8m,+ 8X ) sin r (4u-4y+4iJ. +4X )sin ^ 

+ ( 0 , 0 j k ) ~D ^ g ( 0 , 0 f k ) +Re ( e ika/2 v / 2pdd ( q } l ] k ) ) 

_ . ika/2v^5_dd , n , . . 

- Re(e Djj ( 0, 1 |k) ) 

( 1 . 221 ) 

Now we examine the modes that are two dimensional re- 
presentations. The modes are mixed modes and are not pure- 
ly longitudinal or purely transverse. First consider the 

ZL modes which are a mixture of LA + TO , where TO 
1 z z 

denotes transverse optic polarized in the z direction. 

For these modes 



V°> • b(6 cl +6 a2 ) + c6 c 3 (1.222a) 

V X > - e ika/2 ^[(6 iil+ 6 a2 )-c6 o3 ] . (1.222b) 


Substituting Eqs . (1.222) into Eq . (1.20) and using Eq. 


f 


100 


(1.205) we obtain the secular equation 


I j B i+ B 2 +B 4 +Re [ e 


lka/ayff 


(B 5 +V B 8>]J 


1 

2 


B 1 -f-B 2 -« 4 +Re e 


+ 8 


[ Im ( B 3 _< 


ika/2v / 5’ 
ika/2 V5 b j 


(B k +P c +B q )^ 2 


5 6 8 ' J | 

2 


(1.223) 


In obtaining Eq. (1.223) we have used the facts that 

ika/2v? 


ika/2/? p e ika/2^5 
B l' B 2’ B 4 ,B 5 e ,B 6® 


, and Bg e 
ika/2v/? 


are pure- 
ly real and that B 3 and B.^e“‘*“' *" v " are purely Imaginary. 
The upper sign in Eq. (1.223) is the optical branch and the 
lower sign the acoustic branch. 

Similarly, the ZLj modes are a mixture of LO + TA z . 
For these modes 


V 0) ■ b( 6 al + 6 a 2 ) * c6 *3 

V 1 ' ‘ slka/Z/5 f- b< ^l' t6 a 2 >+<:6 c3 1 • 


(1.224a) 

(1.224b) 


In a similar fashion to Eq. (1.223) we obtain the secular 
equation 


U) 


, 2 - il 


B„ +B n +B .-Re 


± S “l +D 2™4 


e ik»/2V2 (B5+Bg . Bg) 




B^+Bg-B^-Re 


e i^a/2*/^ (B ^ +B ^ +Bg ) 


8 


~ /_ q ika/2V2 \1: 
Iml Bg+B 7 e I 


Jl 


(1.225) 


t 
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The upper sign In Eq. (1.225) is the optical branch and the 
lover sign the acoustic branch. 

At this point we would like to note two additional 
facts. The Brillouin zone boundary along the [110] direc- 

■■ Ow 

tion is the K point and is located at k - —(3/4, 3/4,0) . 

A 

If we continue past the zone boundary to the point 

•^■(1,1,0) , this point is equivalent by symmetry to the X 

point. The reason for this can be seen by using the pro- 
48 

perties 


cu(k + £, j ) - u!(kj ) 


2 ** h* 2 ** 

u> (S«kj) - u> (kj) 


(1.226) 

(1.227) 


together with the reciprocal lattice vector 


6 - ^L(-l f -l,0) 


(1.228) 


and the rotation about the X 2 axis 



(1.229) 


The <j vector takes ^(1,1,0) into 2Tr/a(0,0,l) . The 

fit 

2<jr 

rotation about the xj axis takes —(0,0,1) into 
— (1,0,0) which is the X point. 



I 


t 


t 


fill! DIRECTION 

Along the [HI] direction the wavevector £ is 


-» k „ 

k ■ JS le i 


'2 + *3 ) 


(1.230) 


Direct computation using the model developed in previous 
sections allows one to write the dynamical matrix along 
[111] as 



where and C 2 are real and we have defined 

C 1 - A + Djj(0,0|k) 

C 2 “ B + D^2 (0 f 0 |k) 

c 3 - c + Djj(0,l|k) 

C 4 ■ D + D^(0,l|k) 


(1.231) 

(1.232a) 

(1.232b) 

(1.232c) 

(1.232d) 


and 


A 


1 

M 


4a+8n '+4X '+10c 
+(8n+4\+8u *+4X * )s in 


2 


ka 2c ka 
2v / 3 + ~3~ cos "J3 


(1.233a) 
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B - i(4P+2c/3)»iti 2 


C - - |r(a4Bc/3)(1^3e“ ika/ ^) 


[+(2n'+\') (2+e“ ika/v ^+e" 2ika/v ^ 


D - •g[’(-e+4c/3)(l-©“ ika/v/ ^)+2i/' + (26'-i/')e , ' ika/v/2 


-C 26 '+i/')e“ 2lka/ ‘^ r 


Consider first the A modes which are longitudinal. 
Noting that along [111] 


ik*3(l) -J5 i ka/4 

1 e 


we have 


W ( 0 ) - b ( 5 . +6 _+6 ,) 

a od a2 a3 


. e '/3 - c(6 rf *6 o2 *6 a3 , . 


Substituting Eqs. (1.235) into Eq. (1.20) and using Eq. 
(1.231) we obtain the secular equation 


* _ * 


OJ^* - (C 1 ^C 2 ) ± [ (C 3 +2C 4 ) (C 3 +2C 4 ) ] ? . (1.236) 


At the L point of the Brillouin zone boundary 
£ * ~(£,i,£) and C~ and C. become purely real. We 


then have 
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2 

UJ LAL 


^/6a-4g4«n+Sl/+4X+4u'+2X '+4l/ 
\ + 8u'+4X'+64c/3 


'-86 



+ «ga )+ vga)> 


(1.23 7a) 


and 


where 


2 

“iOL 


i/2a+40+«ti.-^i / +4X+12u '+6X '-4V '+86 ' 

1 

+ 8u'+4X # 


2 

4t (p« +2p«) 22 _22 

cM— < T I?( L )-V?J (L) ) (1.237b) 


T 11 (L) " 2-30147423 
V^(L) - -4.45445881 . 


(1.238a) 

(1.238b) 


Note that the dipoles push the LA mode down and the LO mode 
up. 

Now we proceed to the transverse modes. These modes 
are labeled Ag and are doubly degenerate. One eigen- 
vector is 


w (0) 
a 




) 


V 1 * 


JS 1 ka/4 


b < 6 «l- 6 «2 5 


(1.239a) 


(1.239b) 


Substituting Eqs. (1.239) into Eq. (1.20) and using 
(1.231) we obtain the secular equation 


«•£ ± - ((^-Cj) ± [ (C*-C*)(C 3 -C 4 )P . 
3 


At the L point we have 


and 


where 


°4aL ’ ^/2a-20+8u-4i/+4X+12u'+6X'*2i/' 


-46 *+8y, *+4X *+6c 
2 


4w( Pi -P 2 ) 11 11 

+ iTC (V 12 (L) r i2 (L)) 


a^OL m ^/6a+20+8^-4t'+4\+4u '+2\ '-2v +46 


+8n'+4\ 4 '+34c/3 


_2 I (3P^P 2 ) ( t 22 ( l )+ v 22 ( 


0 


11 v "' ’ll 


L)) 


,11 


T 12 (L) " -3.15523204 


vJg(L) - -6.53319856 
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Eq. 

(1.240) 


(1.241a) 


(1.241b) 

(1.242a) 

(1.242b) 


From Eqs. (1.241) we see that the dipoles push the TA mode 
down and the TO mode up. 



CONDITION OF MINIMUM POTENTIAL ENERGY 


Since we have expanded the potential energy about ita 
minimum value, we must Insure that this condition holds. 

The static potential energy can be written as 

*0 - •y(2N)(4o 1 (r 0 ) + 12© 2 (r 2 > + 12© 3 (r 3 ) + 6© 4 (a)) , 

(1.243) 


where N is the number of unit cells in the crystal. In 
writing Eq. (1.243) we have noted that the dipole terms 
give no contribution to the static potential energy, and 
neither do the angle bending terms. The condition of mini- 
mum potential energy requires 



(1.244) 


Using Eqs. (1 . 9 ) , (1 . 69 ) , (1.81), (1.89) and (1.98) we 
obtain 

a-8 + 8u-8^ + lln'-lly' + 8u' " 0 , (1.245a) 


or 


»i<V 8 «2 (r 2> 1X »3<V 

r 0 r 2 r 3 

M. 


+ 


8© ' (a) 
- 


0 . (1.245b) 


Eqs. (1.245) give a relation between the force constant 
parameters that must be satisfied for our model. Note that 
this condition is trivially satisfied if all first order 
potential derivatives are zero. 
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PHONON DISPERSION CURVES AND PITS TO EXPERIMENTAL DATA FOR 


SILICON 

In this secti'n we present two basic models. First 
we examine a four Darameter model that illustrates the 
effect of the dipoles on the phonon dispersion curvt 
Second we present a fit to experimental data using a ten 
parameter model. 

FOUR PARAMETER MODEL 

This was the first model used when investigating the 
effect of the dipoles on the phonon dispersion curves. In 
this model we take as parameters a m & , \l m v , o and p 1 , 
with all other parameters zero. This corresponds to purely 
central force first and second neighbor interactions, angle 
bending and what we will call diagonal nonlocal dipoles, 
since the matrices Eq. (1.150) are diagonal for this case. 
The experimental data’ fit are C^, C 12 , and Wp AX • 

The parameters determined in this manner are given in 
Table 2. An interesting feature of this model is that the 
short-range force constants are completely determined by 
C ll' C 12 and ' The sin 8 le dipole parameter is de- 
termined by • Since the dipole interaction does not 

affect C 1;L , C 12 f C 44 an d , a comparison of this 

model to one using the same short-range force constants 
but no dipoles can be made. 
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Table 2 - Four Parameter Dipole Model. The experimental 
values of C n , C 12 , and have been used to de- 

termine the parameters. P t - Zje where e is the magni- 
tude of the electron charge in C.G.S. units. 

a • 0 • 3 .672 x 10 4 dyne/cm 

3 

n - v - 3.080 x10 dyne/cm 
a • 7.163 x 10 3 dyne/cm 
- 0.8518 

®'(r Q ) • 1 . 102 x 10^dyne/cm 
® 2 (r 2 ) “ 6 . 161 x 10 3 dyne/cm 

Calculated value of C 44 

11 2 

- 7.476 x10 dyne/cm 
44 

Experimental value of C 44 

11 2 

r. - 7.963 x 10 dyne/cm 
44 
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In Fig. 8 we have plotted th® phonon dispersion curves 
along symmetry directions for the four par am® ter dipole 
model. Also plotted for comparison are the dispersion 
curves for the same short-range force constants but no di- 
poles. Note the dramatic lowering of the TA modes near 
the zone boundaries along the [100] and [110] directions 
due to the dipoles. The TA [111] modes are also lowered, 
but not nearly as much as the [100] modes. The major fea- 
ture of the dipole model is that It allows lowering of the 
TA modes while maintaining a high value of . (Note 

that the calculated value of in this model is within 

6$ of the experimental value.) 

Although the dipoles improve agreement with experi- 
ment, the basic shortcoming of them is that the TAX modes 
are lowered much more than the TAL modes. This is a pro- 
blem with the angular variation of the dipoles which is not 
remedied by addition of the nondiagonal parameter P2 $ 
since the quantity (p 1 -p 2 ) enters into both Eq. (1.200) 
and Eq. (1.241a). It is possible that extending the dipole 
model to second neighbor nonlocality could remedy this de- 
ficiency . 

50 

From the expressions given by Lax we see that his 
quadru poles tend to restore the balance between the TAX 
and TAL lowering. However the quadrupoles also affect 
the elastic constants and Raman frequency. It is the opin- 
ion of this author that including dipole-quadrupole and 
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quadrupole-quadrupole Interact lone would remedy the abort- 
coning of the angular variation of the dlpolee. Thia could 
be done with the addition of only one more parameter for 
the quadrupolea, aa in the work of Lax. The dipole- 
quadrupole Interaction ahould be interacting alnce it la a 
mixture of nonlocal dlpolea with local quadrupolea. The 
difficulty with this approach la application of the Ewald 
method twice more for the long-range auma . 


Ill 


Figure 8. Four parameter dipole model as described 
in the text (solid lines). Circles are experimental data 
(Ref. 49). The dashed lines are the same short-range force 
constants but with no dipoles. 
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TEN PARAMETER MODEL 


In an attempt to Improve the four parameter dipole 
model of the previous section, we have included third and 
fourth neighbor central potential interactions and the non- 
diagonal dipole parameter P 2 . The first derivatives of 
the potentials are also allowed to be nonzero, subject only 
to the condition Eq. (1.245b). The ten independent para- 
meters are ®{(r 0 ) , *£(r 2 ) * 0 3 (r 3^ » ® 4 < a ) » *l (l V * *2 (z V * 
0g(r 3 ) , 0 , P 1 and p 2 . ($ 4 ( a) is determined by Eq. 

(1.245b)). The parameters were determined in the follow- 

2 

ing manner. Using the analytic expressions for Mui^ > 

aC ll ' aC 12 ' MaJ LAX ’ ^AX ' Mui TOX ' ^AL ’ MaJ LAL ’ ^OL ' 

2 

MtOp OL , and aC 44 we performed a nine parameter, weighted, 
nonlinear least square fit to the experimental values of 
these quantities for a given assumed ratio of p 2 //p i • The 
ratio p 2 /P \ was then varied to obtain best agreement with 
the mode near the K point consistent with the smal- 

lest value of least squares of the above quantities. 

The method of solution of the set of nonlinear equa- 
tions is the Newton Raphson iteration with the zero-order 
parameter values determined from a linear least square fit 
to all the preceding experimental quantities except C 44 . 

A typical nonlinear fit converges to an accuracy of ten 
decimals in ten iterations, the value of C 44 changing by 
almost a factor of two from the value computed using the 
zero-order parameters. 
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The entire fitting procedure is made somewhat sub- 
jective by the weighting factors which are necessary to ob- 
tain a reasonable fit to the data. These weights had to be 
Input by hand using the trial and error method and the 
criterion that the percent errors in the frequencies be 
nearly the same. Obtaining such fits to the data is ted- 
ious and time consuming at best. 

In spite of the above mentioned difficulties, we have 
obtained a reasonably accurate fit to the neutron scatter- 
ing data plus elastic constants. The parameters for this 
fit are given in Table 3. In Table 4 we list the potential 
derivatives for this fit. Note the rapid fall off of the 
potential derivatives by the fourth neighbor distance. The 
dispersion curves for the model with these parameter values 
are plotted in Fig. 9. The model reproduces the low-lying 
acoustic modes rather well and the elastic constants are 
within 13$ of the experimental values. Though this is not 
stunning, it is a marked improvement over models employing 
only first thru fifth neighbor interactions of the type 
used in this thesis. The parameters in Table 3 should not 
be taken as the best fit possible with the ten parameter 
model since the author has not exhausted all the regions of 
parameter space in combination with all possible weights in 
the fitting procedure. 



115 


Table 3 - Ten Parameter Dipole Model, “ z t e and 

Pj “ Zje where e is the magnitude of the electron charge 

in C.G.S. units. 

a - 4. 7008 x 10 4 dyne/cm 
5 - 3.0364 xlO 4 

p. ■ 5 . 4455 x 10 3 
v - 4.7303 x 10 3 
X - (^-lO - 7.1520 x 10 2 
a - 4.3265 x 10 3 
u' - -2.2679 x 10 3 
v' - -8.0545 x 10 

# 

X ' - (n'+8i/') - -2.9122 xlO 3 
b' - 31/ ' - -2.4164 x 10 2 
X " - 9.9978 x 10 

u" - - +8n-8i/ +llp # -ll^') - 2.1190 xlO 2 

x 1 - 0.804905 z 2 - -0.201226 z 2 /x \ " “°- 25 


Elastic Constants 



Calculated 

Experimental 

11 

1.4352 x 10* 2 dyne/cm 2 

1.6573 x 10 12 

12 

7.2011 x 10 11 

6.3937 x 10 11 

44 

7.0152 x 10 11 

7.9625 x 10 11 


* (McSkimin, ref. 49) 


Table 4 - Potential Derivatives for the Ten Parameter 


Dipole 

Model. 


■ (a+ 2 e ) - 1 . 0774 x lO^dyne/cm 

®l (r 0 > 

- (a-0) - 1.6644 x 10 4 

r o 

^ 2 (r 2 ) 

- (n+iO - 1.0176 x 10 4 

* 2 (r 2 > 

- X - 7.1520 x 10 2 

r 2 

®3 (r 3> 

- (n'+ 10v') - -3.073 xl0 3 

®3 <r 3 ) 

- (n'- v') - -2.1873 x 10 3 

r 3 

© 4 (a) - 

- 9.9978 x 10 

04 (a) 

a 

n" - 2.1190 x 10 2 


a 
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Figure 9. Ten parameter dipole model as described in 
the text (solid lines). Circles are experiment* 1 data 
(Ref. 49). 
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Chapter 2 


Thermal Expansion of the Diamond Structure 


"For what is a man profited, if he 
shall gain the whole world, and 
lose his own soul? or what shall 
a man give in exchange for his 
soul?" 

Matthew 16:26 


INTRODUCTION 

In this chapter we focus on aspects of the anharmonlc 
lattice dynamics of the diamond structure that are neces- 
sary to calculate the thermal expansion First we review 
some relevant thermodynamic formulas and then develop the 
statistical mechanics treatment of thermal expansion to 
lowest order in the anharmonic ity . Next we present an 
anharmonic model which is a consistent extension of the 
harmonic model of Chapter 1. In the process of developing 
this model we found analytic expressions that explain in a 
simple way how the negative Gruneisen parameters arise. 
Analytic expressions for the mode Gruneisen parameters for 
several modes are presented. A fit to t.ie experimental 
data available for silicon is made and dispersion curves 
for the mode gammas along symmetry directions are presented 
The model is then used to calculate the volume thermal ex- 
pansion coefficient of silicon between 5 and 1700K. Also 
calculated are the thermal strain, zero-point strain and 
zero-point phonon pressure. 
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THERMODYNAMICS REVIEW 


The purpose of this section is to define some thermo- 
dynamic quantities and derive some formulas that are useful 
in connection with thermal expansion. We begin with the 
Helmholtz free energy for the crystal 3, 

F(T,V) - E-TS , (2.1) 

where T is the temperature, V the volume, E the in- 
ternal energy and S the entropy of the crystal. Taking 
the total derivative of Eq. (2.1) we find 

dF - |||j dT + |!^J dV - dE - TdS - SdT . (2.2) 

2 

Using the fact that for a reversible process 

TdS - dE + PdV , (2.3) 


where P is the pressure, we obtain from Eq. (2.2) that 


and 


p - - (w) T • <2 - 4a > 

S - - (f ) y • (2.4b) 


Equation (2.4a) is the equation of state since it relates 
the pressure to the temperature and volume. Note that if 
the pressure is zero, then the crystal volume is determined 
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r 
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from Eq. (2.4a) by setting to zero the derivative of the 
Helmholtz free energy with respect to volume. 

3 

The volume expansion coefficient is defined as 

a ■ ?(lf) p • <2-!S 

Since 4 

(#), - • a,©, • 


we have upon using Eqs. (2.4a) and (2.6) that 


i/dv\ 

01 " V\dP/ T [3T\dV/. 


(2.7) 


Using the definition of the isothermal bulk modulus* 


8 ' * VI $) T > 


( 2 . 8 ) 


we obtain using Eq. (2.4a) that 


B 


- v @) 


(2.9) 


Furthermore, since 



(2 10 ) 
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we have that 



( 2 . 11 ) 


so that Eq. (2.7) becomes 


a - 



( 2 . 12 ) 


Note that the order of differentiation is interchangeable 
in Eq. (2.12). Equations (2.4a^, (2.9) and (2.12) clearly 
show the central role of the Helmholtz free energy. 


I 


STATISTICAL MECHANICS TREATMENT 

The Helmholtz free energy is given in terms of the 

g 

partition function Z by 


F - - k fl T Sm Z 


(2.13) 


where 


7 


Z « Tr e~ eH . (2.14) 

The trace in Eq. (2.14) is over a complete set of states , 
kg is Boltzmann’s constant, 3 - 1/kgT and H is the 
Hamiltonian of the system. For our purpose here, we con- 
sider the Hamiltonian 


H - T + 4 , (2.15) 

with T the kinetic energy given by Eq. (1.15) and t the 
potential energy as a function of nuclear displacements. 
Expanding about the configuration of minimum potential 
energy we have 

t - * o (V o ) +| E t a p(£Kir<')U a (£K)llp(£'K') 

£ 'k'3 

£ K Q 

X'tc'3 

£ 'V'y 

+ • * * . 

(2.16) 
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Here V denotes the volume of the configuration of mini- 
o 

mum potential energy, * Q (V 0 ) the potential energy 

of this configuration, the second order force constants are 
defined by Eq. (1.17) and the third order force constants 
are defined by 




d 3 * 

du (!K)du a U 'K ')du TT r K r 7 
Qt P Y 


(2.17) 


In Eq. (2.16) we have only kept cubic anharmonic terms. 

This is sufficient for calculation of the thermal expansion 
to lowest order in the anharmonic ity . 

In what follows, we will term the "bare" crystal as 
the crystal oscillating about the configuration of minimum 
potential energy, characterized by the position vectors 
Eq. (1.5). As we shall see, this is not the physical cry- 
stal since the presence of anharmonic ity causes a finite 
strain even at the absolute zero of temperature. 

Now consider an isotropic homogeneous deformation of 
the crystal. We define new dynamic displacements v^Uk) 
by the equation 

u a (lK) - eR a (jtK) + v^(iK) , (2.18) 

where e is the deformation parameter or strain. The new 
position vectors of the deformed lattice are 
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£ UK) - (1+€)R (Xk) , (2.19) 

a a 

so that the volume of the deformed crystal Is 

V - (1+€) 3 V . (2.20) 

o 

Substituting Eq. (2.18) into Eq. (2.16) we obtain for the 
potential energy 


* . r 0 ♦ ?> + t a + ? 3 , 


( 2 . 21 ) 


with 


?, 


*0 (V o ) + V p » og U«U'«')H o aK)R p (4'K-) 


X'k 'g 


L 4 a0Y (XK|X / K'|l'K 1, )R a (XK)R e (Jl'K')R Y (l // K") , 

X K Ck 1 


X K p 
X K Y 


(2.22a) 


€ S 4 p (XK |X'*')R ( X K ) Vp (X'k') 
Xk<* aP a e 


X'x'g 


+ V E 4 /vflv (XK|X # K-|X i, *C < ')R /u (XK)R r (X'k')V v (X'k'), 
44 XKq aPY a P Y 


X K p 

X K y 


(2.22b) 
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*2 - i p • gUKU'Ov UOv.U'O 

X Kqj 


X'k'0 


and 


+ § 2 * a3Y Utc|jfc'K'|i # K')R a (lK)v g (l'K')v (X # K') , 

XKa 
X' k'3 

X'k'y (2.22c) 


| £ * g (IK !X'k'|X - K*)v (XK) Vfl (X'Ov(X'k') . 

Ka a p y 

X'k '3 

X'k'y (2 . 22d ) 


We further define the perturbations to the static lattice 
energy and force constants as 


$ 


( 2 ) 1 

m _ 


2 2 • ag (XK|X'K')R a (4K)R g (i'<') , 

X'^g 


(2.23a) 


, (3) 

0 


| 2 * agy (XK |X'k' |X'iC # )R a (XK)R g (X'K')R (X"k") , 

Ji KCfc 


l K p 

X'V' Y 


(2.23b) 


> (1) 

a 


(XK) - • ae (XK|X'tc')R B (X'K') , 


X K p 


(2.23c) 


‘o 2) < iK > -7 2. • ev (iK|i'K'|i*K')R (r K -)R v (i'«') , 

X K 8 


, H » 

X K Y 


(2.23d) 


(XK |x'k') - ^£ # « ag ^(XK |x'k' |X"OR y (X'k " ) . (2 . 23e) 


13 0 

Since the deformation is isotropic, the perturbations to 
the force constants will have the bame symmetry properties 
as the respective force constants of the same order. This 
being the case, we have for the diamond structure that 

t (1) (iK) - 0 (2.24a) 

a 

and 

«< 2) (Xk) - 0 , (2.24b) 

so that 

- 0 . (2.24c) 

Equations (2.24) are shown in Appendix F. Thus we can re- 
write the potential energy as 

i - *0 + i S [* ag UK !x'k') 4- € *^ ) (IK |f'tc')]v a Uic) v U'O 
Si^OC 
X'tc'3 

+ | E « a - Y (XK|X'K'|X < 'K 4 ')V a (XK)v 3 (X'K')V Y (X^ < ') . 
XKqj ' 

X'k'3 

X"k' y (2.25) 

Writing the Hamiltonian as 


H - H q + H 1 , 


(2.26) 
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with 

H 0 * T + W + I , S » oe («|4'«')v o (i'()v 0 Ci'K') , 

l K(X 

I'K'fi (2.27a) 

and 

H x - € 2 ^ 2) +€ 3 *^ 3) +| 2 «^ ) UK|X'K')V a (iK)Vg(i'K') 

L Kot 

+ | Z ft \l\")v (itc)Vg (X'k')vU'k') , 

iKa ap 

iv5 (2.27b) 

we see that H Q is the Hamiltonian of the bare harmonic 

crystal in terms of the displacements v^(Ik) , and H.^ 

contains the anharmonic effects. Note that all the strain 

dependence is contained in the H x term. To lowest order 

8 

in the Helmholtz free energy is given by 

F - F Q + (H^q , (2.28) 


where for any operator M 



(2.29) 

(2.30) 


and 



» 
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-$H, 


Z 0 " Tre 


(2.31) 


Using Eq. (2.27b) we obtain 


<»1>0 - 


2,(2) . ,3,(3) , € 

€ w q + € * 


0 +f s l<J ) U«|l'«') X 

0 2 Xk<* 

x'k 'g 


x <V a (XK)v g U'K ')> 0 , 


(2.32) 


where we have used the fact that' 


<v a (XK) Vj3 (X'tc , )v y (X'k # )) 0 - 0 . (2.33) 


Consider next the normal coordinate transformation 


10 


V a (XK) 


E 

k3 


. e (K|2j)e 111 ' RU) A i * j , 

2NMuj(kj ) 1 “ 3 


(2.34) 


where Ar» is the phonon field operator given in terms of 
the usual creation and destruction operators by 


A 2j ' b itj + b -kj 


(2.35) 


The BZ on top of the sum in Eq. (2.34) is a reminder that 
the sum is restricted to the first Brillouin zone and N 
is the number of unit cells in the crystal. We also have 
the expectation value** 

<A Sj A k'j' ) o ’ 6 k,-S' J jj' (2S 2j +l) ’ (2 ’ 36) 


where 
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(2.37) 


Using Eqs. (2.34) and (2.36) we obtain 


BZ 

(V a (iK )v fj (l'K ')) 0 - L 


-e ft (K |kj)eg (* ' jkj) x 


2NMu)(kJ 

!?• (£(!')-£(*)), „- 


x e 


(2n_^ +1) . 
kj 


(2.38) 


Substituting Eq . (2.38) into Eq. (2.32) yields 


<«1>0 - *\ 2) + e3 *o 3) 


+ ^f £ * 


4NM 


kJ l' K'& 


^ e (< |kj )e a (K ' | kj ) x 
iKo aj(kj) “ p 




(2.39) 


Now we note that 


2 .UK E^>0CK'|S> , 

Jl it 


(2.40) 


where we have defined the perturbation to the dynamical 
matrix in analogy with Eq. (1.25), 

d£* ) (K , c'|£) -is •^J ) (0K |nOe i * , * <m) . (2.41) 


M “ ap 
m 
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Thus Eq. (2.39) becomes 


<%>0 - 2) ^ *^ 3) 


. BZ (2r.+ .+1) j . /« . . 

+ &Z 2 e*<« |lfj)D5.| ) (KK'|S)e a CK' |Cj ) . 

4 *, «a a(Jj) a " P 

J (2.42) 


Defining 


- £ e*(tc!?j)D^ ) (KK'|^)e g (K'|itJ) * 


(2.43) 


Ka 

k'P 


we have In more compact form 


< K 1>0 


.2,(2) 3 (3) eft B? "l (k J ) - ,i, 

0 0 2 


(2.44) 


Since H Q does not depend on the strain, to lowest order 
in we have from Eq. (2.28) that 


dF m a<H i^o 
de de 


(2.45) 


so that 


dF 

de 


¥• - 2et 


(2) , « 2. (3 ) . n r 
0 + J€ *0 + 1 ^ 


n I? “i< k J> - 


Sj a ' ,(kj) 


(njj+i) . (2.46) 


To proceed to volume derivatives we use Eq. (2.20) so that 
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3V 0 (l+€)‘ 


dF 

■*€ ‘ 


(2.47) 


Using Eq. (2.4a) we obtain 


3 V n (l+€ ) L 


2e*l 2) -f-3€ 2 * 


) J J 


(2.48) 


From Eq. (2.9) the bulk modulus of the bare crystal in the 
harmonic approximation is given by 


B o ■ ~5VJ“ ' 


(2.49) 


so that to lowest order in € the equation of state is 


P - 3B 0 £ 


u£(kj) 


■§ Y~ u — “ — 

o ^ w(kj) 


(n icj + 


(2.50) 


Defining the mode Griineisen parameter by 


Y (kj ) 


u^(itj) 

6a) 2 (kj ) 


(2.51) 


we have 


P - - 3Bq£ + y L oj(kj )v (£j) (nj^j + £) . 
0 it j 


(2.52) 


From Eq. (2.52) we see that even at the absolute zero of 
temperature and zero pressure there is a zero- point strain 
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given by 


BZ -# 

£ 0 • 3B^f 5a ^ a Y(f J ) • (2 ' 53) 

“ J 

This can be thought of in terns of the pressure required 
to make the strain zero at T - 0 , 


p 0 - 4 - ? ba&LL y ( 5 j ) . ( 2 . 54 ) 

P Q is thus a zero- point phonon pressure. 

Proceeding now to the thermal expansion coefficient, 
we have from Eqs. (2.12) and (2.47) that 


a 


1 

3BV q (1+€) 2 



» 


c 


(2.55) 


or using Eq. (2.46) 


ft BZ u£(kj) d t , v 

O' 3 - 6 ^ 3T n k1 

6BV Q (l +€ ) A gj w(kj) kJ 


(2.56) 


Now 


+ *> ’ ^^| Sinh 
B 


2 


I- 


J| 


(2.57) 


Defining the Einstein specific heat function as 


'^(a:) • j sinh t Pftaj/2 ] | 


(2.58) 


Eq. (2.56) becomes with the aid of Eq. (2.51) to lowest 
order in e 


i BZ 

pTT” £ Y (llj )C p (u>(icj ) ) . 
E 


(2.59) 


Eq. (2.59) is the same as the result of the "quasi-harmonid' 
12 

approximation, at least to lowest order in e . This 
approximation is useful to find the meaning of the mode 
GrUneisen parameters. In the quasi-harmonic approximation, 
the phonon frequencies are given to lowest order in e by 


2 2 (kj ) - u ) 2 (itj ) + eu£ (k j ) + • • • , 


(2.60) 


with ui^(kj) defined by Eq. (2.43). To make connection 
with experiment, one defines the mode gamma by 


Y(iTj) - - 

d QnV 


(2.60) 


or 


Y (kj ) - - 


V 

iiFcirjT 


dS 2 (gl) 

dV 


(2.61) 


Using Eqs. (2.20) and (2.60) we obtain 


Y (kj ) - - J (ltj ) 

6 w (kj) 


(2.62) 


so that to lowest order in € 
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- «f(kj> 

Y(kJ) - - -V— . (2.63) 

6ur (k j ) 

Equation (2.63) agrees with the definition Eq. (2.51). To 
proceed one atep further, experiment* measure the pressure 
dependence of phonon frequencies for some modes and thus 
one can obtain experimental information about the mode 

gammas . 



ANHARMONIC MODEL 


In order to calculate the thermal expansion using Eq. 
(2.59), it is necessary to have a model for the mode gammas. 
The earliest model of GrUnelsen was to assume that y(kj) 
was a constant, the same for each mode, and thus the ther- 
mal expansion would have essentially the same temperature 
dependence as the specific heat. Th^s is known as 
Griineisen's rule. It is a reasonable approximation for 
some materials, but since the specific heat and bulk modu- 
lus are always positive, it is completely inadequate to 
describe materials with negative thermal expansion that 
have a change in sign of the expansion coefficient. 

It is interesting to note that a linear monatomic chain 

of atoms with nearest neighbor interactions gives a mode 

14 

gamma that is the same for each mode. Thus we have one 
example where Griineisen's rule is rigorously justified. 
However, the inclusion of second neighbor interactions in- 
troduces dispersion so that the mode gammas become wave- 
length dependent even in this simple, one-dimensional 

. i 15 

model . 

Several discussions of negative thermal expansion have 
appeared in the literature. Experimentally, negative ther- 
mal expansion has been observed in the tetrahedrally bonded 
solids 16 C, Si, Ge, GaAs, GaSb, CdSe, CdS, CdTe, Agl, 

InSb, CuCl, CuF, Cul, HgTe, ZnO, ZnSe, ZnS, ZnTe, hexagonal 
ice and CuInTe 2 • Negative thermal expansion has also 
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been observed in the tetrahedral glasses Si0 2 , GeOg , 

BeF 2 and Si0 2 + x^Na 2 0 for x less than 0.2 . Nega- 
tive thermal expansion also occurs In other crystal struc- 
tures (fee RbBr and Rbl , for example) but is relative- 
ly rare in comparison with the tetrahedral solids. 

On the theoretical side, the occurence of negative 
thermal expansion may be understood from Eq. (2.59) if the 
mode Griineisen parameters are neg&tix'e for enough modes 
that contribute a large weight to the sum in Eq. (2.59). 
Although several discussions of negative Griineisen para- 
meters have appeared in the literature, none of them have 
been very realistic or explicit in explaining how the ne- 
gative mode gammas arise. For this reason we ise a model 
with sufficient generality to explain the origin of the 
negative Griineisen parameters and present analytic expres- 
sions for them for several modes. 

The anharmonic model we use is a consistent extension 

of the harmonic model developed in Chapter 1. The reason 

for mentioning this is because previous calculations for 
17 

silicon have used a shell model for the harmonic proper- 
ties and a rigid ion model for the anharmonic properties, 
with no relation between the two. In fact, the first and 
second order potential derivatives which enter the harmonic 
model also strongly affect the anharmonic model, as we 


shall see. 
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In order to compute the mode Griineisen parameters, 
once we have the harmonic and anharmonlc models, we only 
need compute the perturbation to the dynamical matrix Eqs. 
(2.41) and (2.23e). Since the perturbat Ions to the force 
constants have the same symmetry properties as the second 
order force constants, once we have obtained the form of 
the dynamical matrix for the harmonic model, the form of 
the perturbation to the dynamical matrix can be written 
down Immediately. The actual values of the perturbations 
must be calculated from physical considerations. 


CENTRAL POTENTIAL CONTRIBUTIONS 


For two body central potential interactions we can 
write tbe cubic enharmonic contribution to the potential 
energy as 


t 


3C 


A Z Z Z <t - (ikU'k')u UkIi'k') x 

12 IK i'K- oBy “ 


x U g (IkJx'ic')u (JlK|l # tc') . 


(2.64) 


This follows from Eqs. (1.54) and (1.56) and the coeffici- 
ents (IK | i'K ' ) are defined by Eq. (1.57c). By per- 

forming an isotropic, homogeneous deformation Eq . (2.18), 
one can show that the perturbations to the force constants 
0 ^Uk|x'k') defined by Eq. (1.57b) are 

®^ } Uk|X\'') - Z 0 a&y (IK U'k')R y UK jX'O . (2.65) 


FIRST NEIGHBOR CONTRIBUTION 

By direct computation using Eq. (2.65), we have in 
complete analogy with Eq. (1.68) that 


® (1) (0,0|6 lf l) 



( 1 ) 

p(l) 

B (1> \ 

( 1 ) 

* (1) 


( 1 ) 

e (1) 

•«> 
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0 (1) (O,O|6 2 ,1) 


*d> 

*<*> a ( x > -p (1) 

- P (D ■*<*> 


® ( 1 ) ( 0 , 0|6 , 1 ) - 


* (1) -a (1) -a (1)> 

-3 (1 > * (1 > p (1) 

. 0 U) p (l) a (l), 


t ( 1 ) ( 0 , 0 | 6 4 , 1 ) 


a (1) -e (1) 0 (1) ' 

■p^) a (1) -a (1) 

0 (1) -g (1) a (l> 


where we have defined 


a 


( 1 ) 


O m 


2»i< r o> 


T»I ( r 0 ) + f»l< r o)-3 r 


(X) 


°»i ‘V -J«i (r o ) + 3 F~ 4 i <r o ) 


Note also that 


<a (1 >-S (1 >> 


D[ (r ) 

, # / \ o 

»l (r o ) -— — 


Using Eq. (2.41) and the analogy with Eq. (1.25) 
write down the first neighbor, central potential 
tion to the dynamical matrix 


( 2 . 66 ) 


(2.67a) 

(2.67b) 


( 2 . 68 ) 

we can 
perturba- 
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d‘J u (0,l|k) - 




' 5 p. • 1 , 


1 0 o' 


s(l)l 


(0,°!^) - ^2LllJo 10. 


Vo 0 1 


The are defined by Eq . (1.62). 


SECOND NEIGHBOR CONTRIBUTION 


3> <1) (o,o|^ 1 ,o> - 


® (1) (0,0|n 2 ,0) - 


« (1) (0,0 |^„, 0 ) - 


( 1 ) 


( 1 ) 


( 1 ) 


(1.80) 

are 


0 ) 

, <l> 

0 

0 

x (1 ) 

0 

0 ^ 

, (1) 



(1) 

0 

„ (1 > 

x (1) 

0 

0 

, (l) 


® ( 1 ) ( 0 , 0 |u 4 , 0 ) - 


*<!> -„ (1) 0 

-, (1) ^ (1) 0 


0 


( 1 ) 


9 ( 1 ) ( 0 , 0 |^ 5 , 0 ) 


® ( 1 ) ( 0 , 0 |^ 6 , 0 ) - 


(1) 

0 

0 1 


u (1) 

- (1) 


-, (1) 

. a) i 

(1) 

0 



x (1 > 

0 

(1) 

0 

u (1) 


( 2 . 70 ) 
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Here we have defined 

" 2^ r 2®2^ r 2^ + ®2^ r 2^ " ^2^ r 2^ r 2^ 
l/(1> “ | (r 2 0 2 (r 2 ) ~*2 (r 2 ) + ^2 (r 2 )/r 2 ) (2 

\ (1) - (^ (1) -!/ (1) ) - (0'(r 2 ) -02 (r 2 )/r 2 } * (2 

Thus we find that 

D ^ )2 ° (0 ,i|k ) - 0 (2 

and 

o 6 

D^ )2 (0,0|lc)-f 2 0^ ) (O,O|n 1 ,O)[l -cosCi?.^)] 

(2 


The n i 


71a) 

71b) 

71c) 

72a) 

72b) 


are defined by Eq. (1.75). 


THIRD NEIGHBOR CONTRIBUTION 


The matrices analogous to Eq. (1.88) are 


£ (1) (0,0 | T x , 1 ) 


$ ( 1 ) (o,o|t 2 ,i> 


® (1) (0,0|t 3 ,1) 


? (1) (0,0|t 4 ,1) 


® (1) (0,0|t 5 ,1) 


0 (1) (O,O|t 6 ,1) 
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£ (1) (0,0jT 7 ,l) 


0 (1) (O f O|T g ,l) 


0 (1 ) (0,0|t 9 ,1) 


® (1) (0,0|t 10 ,1) 


0 (1) (0,0|t 11# 1) 


0 (1) (0,0|t 12 ,1) 



(2.73) 


Here we have defined 
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“' C1) - TTt r 3«3 > *'3 > °* 3 ( r 3 > - 1 °®3 < r 3 > /r 3 1 

v’ ll> - rc(V? r 3 ) -*3 (r 3 )4 *3 <r 3 )/r 3 1 

6' (1) . 3l/' (1) " n;t 1 ‘3»3 <r 3 ) '®3 <r 3 )+ ®3 <r 3 )/r 3 


(1) . (n' ( 1 ) + 8 l'' a) ) - 


2 * 3 * p 3 * 


9 «/. v 2 * " \ _ — 

ll r 3®3 (r 3 ) H®3 (r 3 ) 11 r. 


Thus we find that 

,o 


D^ )3 (0, 1 |k) - 


,12 


and 


D (1)3 


(0,0|k) - 


(8u 


(1) 


+4X 

M 


( 1 ) 


/I 0 0 \ 

4 o i o ) 
\ 0 0 if 


The ^ are defined by Eq. (1.85). 


(2.74a) 

(2.74b) 

) 

(2.74c) 
(2 . 74d) 

-+ 

T i 

(2.75a) 

(2.76b) 



r 


FOURTH NEIGHBOR CONTRIBUTION 

The matrices analogous to Eq. (1.97) are 

0 (1) (O,O|X 1 ,O) -/x' (1) 0 ° \ 

U (1) o ) 

\o o 

$ ( 1 ) (o,o|\ 2 ,o) - L' (1) 0 0 \ 

(o x'< X) 0 ) 

\o o u' (1) / 

0 (1) (O,O|X„,O) - /n #(1) 0 0 \ 

o ,*(» o 

\o 0 x' (1) /. (2.76) 


Here we have def lned 


\ ,(l) - a«'<a) 

(2.77a) 

a' (1) - (»<(a)-i *;(*)] . 

(2.77b) 

Thus we find that 


D ^e ,4 ° (0 > 1 l J) ' 0 > 

(2.78) 


and 
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1S1 


D ^ )4 < 0 , 0 |£) 


,( 1 ) 4 * 


ora 


<0,0|if)6 ae 


(2.79) 


with 


D ii )4 ° (o »°i^ ) 

D 22 )4 ° (0 » 0 I^ ) 

D 33 )4 ° (0 *°^ ) 



Here „ kg , k^ are the cartesian components of the wave- 
vector defined by Eq. (1.101). 



ANGLE BENDING CONTRIBUTION 


As was mentioned in Chaptsr 1 , ths angle bending 
potential energy Eq. (1.104) contains enharmonic terms 
when expressed in terms of the displacements u (IK) . In 
order to calculate the perturbation to the dynamical matrix, 
Eq. (2.41), caused by the angle bending anharmonlcity , it is 
only necessary to obtain the dynamical matrix with the de- 
formation Eq. (2.18) present, and then separate out the 
term linear in € . This approach is chosen since the 
angle bending contribution to the dynamical matrix has al- 
ready been obtained in the strict harmonic approximation. 

Using Eqs . (1.114) and (2.18) we have 

x(iK|i'K') - (1+€)$UK ]i'*') + vUK | f'O . (2.81) 

Since the deformation is isotropic, it is clear from Ec»s. 
(1.113) and (1.115) that the angles when all the v(Xtc) 
are zero are the same as the equilibrium angles. By using 
Eq. (2.81) and carrying out an expansion similar to the one 
leading to Eq. (1.117), we obtain the change in angle cor- 
rect to first order in the displacements v(l k) and all 
orders in € 

1 a 

(2.82) 
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Here the r\ a (L*.\ L* k' \ l' *.' ) are defined by Eq. (1.118) in 

term* of the undeformed or "bare" position vectors lt(4K) . 

Equation (2.82) differs from Eq. (1.117) by the factor 

(1+e) and the fact that we are using the displacements 

v(1k) . Since we square the angle changes in Eq. (1.104), 

the only difference between the angle bending potential 

energy which is quadratic in the displacements u(1k) and 

that which is quadratic in the displacements v(4k) is 

-2 

that a factor of (1+c) appears in front of the latter. 
Thus the angle bending contribution to the dynamical matrix 
with the j train present is the same as without the strrin, 
Eqs. (1.124) and (1.125), but with the fcrce constant a 
replaced by a where 

c - 2__ . (2.83) 

(l+€) 2 

The contribution to first order in e is 

c - a + 60^ + **• , (2.84) 

where from Eq. (2.83) we have 

c (1) - - 2a . (2.85) 

Thus the angle bending perturbation to the dynamical 
matrix \kk ' (k) is given by Eqs. (1.124) and 

(1.125) with a replaced by everywhere. Rather 

than reproduce the lengthy expressions which are obtained 
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trivially from Eqs. (1.124) and (1.125), we give an example 
for one element 


.( 1 )( 6 ) 

*12 


4o 


( 1)1 ilc* iic'tg ^3 


ik 


(0,1 k) - - 5¥ - \l->-e -e -e 


-e 


•M. 


( 2 . 86 ) 


Eq. (2.86) is the analog of Eq. (1.125b). The t are 
given by Eq. (1.62a). 

The result Eq. (2.85) is rather interesting because 
it predicts that the angle bending force constant pertur- 
bation is given in terms of the parameter a that is de- 
termined by the harmonic model. Thus this type of anhar- 
monicity requires no additional parameters that cannot be 
determined from the harmonic model. This is in contrast 
to the central potential anharmonicity which requires the 
third derivatives of the potential as additional parameters 
that are not determined by the harmonic model. To our 
knowledge, the result Eq. (2.85) has not been obtained 
previously . 


NONLOCAL DIPOLE CONTRIBUTION 


In order tq compute the 'le contribution to the 
inharmonic ity, we must generalize Eq. (1.129) to Include 
terms quadratic in the displacements so that 


p c (ilt) ■ P <rf U | '!*'0< , 9 u'0 

L K p 


+ * P a g Y (^ K U /K 'l^ #K ')UgU'K')u U'k') . 

1 K 3 


A k y 


(2.87) 


Here we have defined 


P 


q3 


(AK |A'k') 


dp (IK) 

*r*T rri 


9 


o 


(2.88a) 


and 




a 2 p a (iK ) 

durri^Hiu tftpi 

P Y 


(2.88b) 


The total dipole moment of the crystal Eq. (1.130) can then 
be written as 


P 

a 


is £ s P a <AK|A' 
2 IK i' K'j3 A'k\ a0Y 


l'K*)u 0 (l'K')U Y (i"K") , 

(2.89) 


where we have taken account of the fact that the term 
linear in the displacements was shown to vanish in Chapter 
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1 (Eq. (1.140)). Defining 

- 2 P aey (^|i'K'|i'K') , (2.90) 

X K 


we have 


P 


Of 


iS 2 M a (XK | £ K )Uq(XK)u (X K ) 
1 XK 6 X'k'y a0Y 0 Y 


(2.91) 


Eq. (2.91) leads to the two-phonon infra-red spectrum of 

18 

diamond structure crystals. 

If we substitute Eq. (2.87) into the dipole-dipole 
interaction energy, Eqs . (1.153) and (1.1551 the resulting 
expression will contain cubic and quart ic anharmonic terms 
in addition to the harmonic term treated in Chapter 1. The 
perturbation to the dynamical matrix can be determined by 
computing the dynamical matrix with the strain present and 
then separating out the term linear in € . 

Substituting Eq. (2.18) into Eq. (2.87) we obtain 

p a (XK) - €P^ 1 } (XK) + € 2 p£ 2 ) (XK) + 

+ S p o (Xk | X k )v„ (X KX) j - 
x'ic'e ^ p 

+ € 2 p^ } (XK |x'k')v p(X'k') + 

X'k' 0 0,3 p 

+i 2 p (Xk|X'k']x'k')v b (X'k')v (X'k') , 

* l'K'& aPY P Y 

A 11 11 

X K Y 


(2.92) 
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where we have defined 


pi 1 ) (lK) - Z p (Xk|X'k')R b (1'0 , 
1 k'P p p 


(2.93a) 


p' 2) (iK, 


-§ _ E, a P [ , ev (iKU''<'|i'«' ) B p (l' | t')R v (i'K') , 


X K p 
X'k'y 


(2.93b) 


and 


•>ae ,(iK 


X K ) 


X*K 


.W*" 


(i'O. (2.93c) 


At this point we impose infinitesimal translation in- 
variance on the moments p a (XK) . Note that the results 
of Chapter 1, Eqs. (1.149), satisfy infinitesimal transla- 
tion invariance without ever imposing this condition on 
the moments p (XK) . Imposing this requirement on Eq. 

Of 

(2.87) we obtain the conditions 

Z p „(X.|X'k') - 0 , (2.94a) 

X'K' ' 


and 


Z p (Xk I x'k ' | x'k") - . Z p « (xk |x'k'|x"k") - o . 

1 i ' k ' 

(2.94b) 


X K 


Infinitesimal translation invariance of the total dipole 
moment Eq. (1.130) only yields the less restrictive condi- 


tions that 
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Z £ p a <£K|£'K # ) - 0 , (2.95a) 

IK I'K' °® 


and 


Z Z 

XK t'K' 


P^^AKfi'K'll-K') 


E £ 

XK X' K 




(xk |x'k' |x'k') - o . 

(2.95b) 


Using Eqs. (2.94), (1.142) and the transformation pro 
perty under the space group operation { S | v(S)+i$(m) } that 


Pc6Y aK|L ' K 'l L ' K ' ) ■ S . S «uV S YX , V>'X <1K ! r, ''l 1 '*' ) ’ 

\jtlr \ 

(2.96) 


it is easy to show that p^^(XK) and p^ 2 ^(£K) transform 
in the same manner as the first order atomic force con- 
stants. Thus by the same arguments as presented in Appen- 
dix F, p^^CXK) and p £ 2 ^(Xk) are both zero in the dia- 
mond structure. Similarly, one can show using Eqs. (2.94b) 
and (2.96) that the coefficients p^j^ ( Xk | £ ' k ' ) have the 
same symmetry properties as the p^ (XK | X 'k * ) . Thus we 
can use Eq. (1.150) to write 
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5 ( 1 ) ( 0 , 0 | 6 1 , 1 ) 


; (1) (0,0|5 2 ,1) 


p (1) (0,0|6 3 ,l) 


p (1) (0,0|6 4 ,l) 


„(1) 

(1) 


Pi 

p 2 

p 2 

n (1 ) 

B (l) 

„(1) 

p 2 

P 1 

p 2 

_ (l ) 

(1) 

(1) 

p 2 

p 2 

p l 

_ (l ) 

(1) 

(1) 

p i 

p 2 

- p 2 

n (i) 

„(1) 

(1) 

P 2 

p l 

” p 2 

(1) 

(1) 

(1) 

P 2 

- p 2 

p l 

(1) 

( 1 ) 

(1) 

p l 

~ p 2 

-p 2 

(1) 

(1) 

_ ( 1 ) 

p 2 

p l 

p 2 

_ ( 1 ) 

(1) 

(1) 

P 2 

P 2 

p l 

_ ( 1 ) 

(1) 

(1) 

P 1 

- p 2 

p 2 

( 1 ) 

(1) 

(1) 

p 2 

p l 

" p 2 

(1) 

(1) 

(1) 

P 2 

’ p 2 

p l 


. (2.97) 


We also have for the case of nearest neighbor nonlocality 
that 

Pa0 J (0> 0 | 0, 0) - - 4p^ 1) 6 ag . (2.98) 

Here p^ 1 ^ and P 2 *^ are P a r an, eters to be determined and 
have the units of charge. From the preceding remarks and 
Eqs. (2.92), (2.97), (2.98) and (1.49) we can write 
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P a U,0) 


- £ £ P a g(0,0|6 i ,l)(v p U+6 i ,l)-v g (i,0)) 

@ 1*1 

+ | £, P agY U,0|X'K'|X # K # )v g (i'K')v y (i # K') , 

i K $ 

(2.99a) 


X'k'y 


p^' 1 * 


£ £ p vfl (0,0|6.,l)(v fl (l-6.,0)-v g(l f l)) 

0 i-1 ** 1 Pi P 

+ 3 P^Y**' 1 ^* ^ i#K * )v 3 (A ' 

X K p 


3 

X'k'y 


(2.99b) 


where we have defined 


? ag (XK |X' K ') 


( 1 ) 


P a( 3 (XK|X'K') + € p^'UK|X'K') . 


( 2 . 100 ) 


In order to compute the dipole contribution to the 
dynamical matrix with the deformation present, it is only 
necessary to keep the terms linear in the displacements 
v^(XK) in Eqs. (2.99). This is so since the quadratic 
terms give rise to terms in the potential energy that are 
cubic and quart ic in the v q (Xk) . Thus we see that it is 
only necessary to replace p^ and p 2 in Eqs. (1.165) by 
their renormalized counterparts and p 2 where 


( 1 ) 

P x - P x + €Pi , 


(2.101a) 


P 2 * p 2 + €P 2 


( 1 ) 


(2.101b) 
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Separating off the terms linear in e we obtain the per- 
turbation to the dipole contribution to the dynamical mat- 
rix 


,(1 )dd 


(0,0'k) - 




€ Mu 

a o 


T ;i^ k > + 


D^ 1 )dd (o,i!k) 


4w(P i P< 1) +p 2 p{ 1) ) 12 _ 

+ — Hr, — — T ol ( k) 

s o 

8ir P 2 Po 22 -+ 

+ T““(k) , (2.102a) 

£ Mfi 

s o 

8irp i p i 1) ai - 

" “TBT — ae (k) + 

s o 

4*( Pl P^ 1) +P 2 P^ 1)) >JL 2 /!t>% 

8,fr P2P2 22 -*• 

+ - & aB ( k > • (2.102b) 

s o 

The matrices T**(k) . T*g(k) , T^(k) , v“(k) , V^(k) and 

are defined by Eqs . (1.166) and (1.167). In Eqs . 

(2.102) £ is the static dielectric constant, not to be 

s ’ 


*5* 


confused with the strain parameter 


Note the interest- 


ing way in which the harmonic and anharmonic parameters 
both play a role in Eqs. (2.102). 


I 


grPneisen parameters 

In order to compute the Griinelsen parameters, there 
are two alternative but equivalent methods. One method is 
to simply use Eqs. (2.43) and (2.51) and directly compute. 
This method is well suited for numerical work at an arbi- 
trary k point once we have obtained the eigenvectors. 

The other method is to compute the quaslharraonic frequency 

N 

u>(kj ) defined by the eigenvalue equation 


£ 2 (kj)e a (K |£j) - 5 ff a0 (KK'!£)eg(K'|£j) , (2.103) 

tc 3 


with 


0 ag (KK'[k) -D ag (KK # |ii) +€D^ ) (KK'|k) . (2.104) 


It is easy to show using first-order perturbation theory 
that 


d ^ 2 1 ih ) 


Ti 




G-O 


- OjJ(kJ) 


Z e*(K|kj)D^g ) (KK'|k)e 8 (K'|kj) . 

K at p p 

k'3 (2.105) 


Using Eq. (2.51) we can write 


Y (k j ) 


1 

a<2 2 tkjj 

6co 2 (kj ) 

de 


(2.106) 


Thus we see that if we have an analytic expression for the 
quasiharmonic frequency of a given mode, one can 
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differentiate analytically and use Eq. (2.106) to obtain 
the mode gamma . 

The quasiharmonic expressions for the frequencies can 
be easily obtained from the harmonic ones by simply re- 
placing the harmonic model force constants by the quasi- 
harmonic force constants according to the prescription 

- l ae (iK|l'K # ) + . 

(2.107) 

Since we have already obtained expressions for the harmonic 
frequencies of several modes in Chapter 1, we will adopt 
this second method to obtain analytic expressions for the 
mode gammas. Note that both methods give identical re- 
sults but are different computational schemes. 
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flOM DIRECTION 

In this section we list some results that can be 
easily obtained from Chapter 1. As an example, consider 
first the Raman mode. Using Eq. (1.190) we can write 

‘Sra “ a + ieiT' + 8D , (2.108) 


where 


( 1 ) 

or - or + ear , 


. ^ *( 1 ) 
or ■ u + £W , 


r* - x' + €\' a) 


(2.109a) 

(2.109b) 

(2.109c) 


Thus 


- i(8a (1) +-~ c (1) f 16^' (1) +8\' (1) ) , (2.110) 


and using Eq. (2.106) we obtain 


RA 


V(8« (1) 4 a (1) + 16.' (1> +8X' (l) ) 


6M “*A 


( 2 . 111 ) 


This can be written in terms of the potential derivatives 

as 
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RA 



o 




8 »l (r o> 


3i Si 


+4r 3*3 (r 3 ) ‘ rt «3 <r 3 > Sr - *^ 


’) ' 


( 2 . 112 ) 


where we have used Eqs . (2.67a), (2.74) and (2.85). 
At the zoDe boundary wt obtain from Eq. (1.192) 


Y LOX " V LAX 


or 


1 

— 5 — 
3m “lax 


a. a > + ^o‘i> + «» < i) + 4u' <1 ) +a i'<iS 

4,r,, i> , i 1) ,ii 


£ . r o 


T ii< x > 


(2.113) 


LAX 


r o x 2 ^ , 2 ®l (r o ) 

T W + 3 W -3 “ 7 — 


3M “LAX L 

2 $ I (r 0 ) 

+ 2r- 2 * 2 (r 2 ) + 2»'(r 2 ) 


2» '(r,) 

+ r 3*3< r 3> +29 3<r 3 >---?r^- 


( 1 \ 

20 2,rp l P l T ll rx) 

3 c + n T 11 (X) 
8 O 


(2.114) 


From Eq . (1.200) we obtain 


» 
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2a (1) -20 (l > + 44 (1 > + 4X (1) + 4u' (1) ^' (1) 




- 2*** x '+46 #l *'+«0 


e . /”^Pl“* > 2^ |:> l ^” p 2 ^ t 22 VA ’^' fV 23 ^ x ) )j * 


(2.115) 


or In terms of potential derivatives 




3Mw TAX 


r + r 2®2 + 3 ®2 


3 ®2 (r 2 ) . 16 ./«, x . 17 x 

“77 + TT r 3^3 r 3 + TT ®3 (r 3 ) 


17 ®3 (r 3 * 


- 6c 


+ 7 %-(p 1 -p 2 )(p‘ 1 , -p< 1 ) ht“(x)+^(x)) , 

_ SO _ 


(2.116) 


where from Eq. (1.202) we have 


T 22 (X) + V 23 (X) “ “16.28143 . 


(2.117) 


Similarly from Eq. (1.201) we obtain 


3MaJ TOX 


2«< 1 W 1 W 1 W 1 W (1 W (1) 


+ 2,' (1) -46'< 1 >+f a (1) 


+ I^ (P 1 +P 2 )(P 1 1)+P 2 1))(T 22 (X) - V 23 (X) ]J * 

(2.118) 
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In terms of the potential derivatives 


TOX 


2 

3M<4 rox 


* , 3 «2 (r 2> 

+ r 8 » 2 (r 2 )+3» 2 (r 2 ) 

6 /// 27 0 27$ (r« ) 

+ IT r 3®3 (r 3 ) + TT ®3 (r 3> H 7~ ■ 


“ ^ + ^?r(Pl + P2 )( Pl 1)+ P2 1))(T 22 (X) 

s o 


C 


where from Eq. (3.202) we have 


(2.119) 


T 22 (X) “ V 23 (X) " 10.76338 . 


( 2 . 120 ) 


In the elastic region we have 


v (el) 
y LA[100] 


6C 


11 


de. 


ii 


d€ 


Je=*0 


(2.121) 


where from Eq. (1.184) 


a^ lx - (a+4o+8(T+2(r'+9r'+8r / ' ) , (2.122) 


and in addition to Eqs. (2.84) and (2.109) 


m p. + €p. 


( 1 ) 


r - x" + € \' (1) 


(2.123a) 

(2.123b) 
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Using Eqs. (2.121) and (2.122) we obtain 


v (el) 
y LA[100] 


6aC 


1 (a (1) + 4o (1 W (1) + 2,' (1) + 9X' (1) 48r (1) ) , 


11 


(2.124) 


and In terms of potential derivatives 


v (el) 
y LA[100] 


Similarly 


6aC 


11 


O . /// 


2 ^ u 


2 »i ( V 


T ®l (r o> + f »l (r o ) ’3 . 

o 

4C«( r„) 

+ 4r 2*2 (r 2 )+44 2 (r 2) 

, 83 _ x , 38 x 38 ®3 (r 3 ) 

+ IT r 3®3 (r 3 + TT ®3 (r 3 ) "IT —— 


8a® 1(a) - 8cr 


(2.125) 


v (el) 
Y TA[100] 


1 

r^4 

6C 44 

d<£ 


Jc-0 


(2.126) 


where from Eq. (1.199) 


»e 44 - 3+| ? + c-4C-r'+ior'-«c 


(| ?-r+2S'+3p-> 2 

(3 + 1 3+23" +rT 


(2.127) 


and 


r - x + ,x (1 > 


(2.128a) 


169 



- u' + €u' (1) 

(2.128b) 

% 

-» + e0 (1) 

(2.128c) 

V 

- 6' ♦ C6' (1) 

(2 . 128d) 

V 

- + «.*<» . 

(2 . 128e ) 


Differentiating Eq. (2.127) we obtain 


Y 


(el) 

TA[ 100 ] 


6aC 


44 


a a) + | o a) + 4 ( 1 (i> + 4 X (i> +x .<i> +10ll . 


-|(2A (1) -|b (1 >) 


( 1 ) 


(2.129) 


where we have defined 


- 

(| o-S+26 '+3l/') 

(2.130a) 

as 

(a + ^ 0+ 2u +X } 

(2.130b) 

CD . 

(i a< l) -0 (1, +26'< l W (1 >) 

(2.130c) 

d) _ 

(a (D + | o ( 1 W ( 1 ) +X' (1) ) . 

(2 . 130d) 


One can rewrite Eq . (2.129) in terms of potential deriva- 


tives as 
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(el) 

y TA[100] 


6aC 


44 


o ^ */ 


X 0 l (r o ) + f •l (p ~ ) 


2 *i (r o } 


o' 3 r 


+ 2r 2 ®2(r 2 )+64j(r 2 ) f-2- 

la 2«'(r 7 ) 

+ irv; (r 3 , ^3 ,(r 3 ) -V" 

+ 8«;<a) 1 1 a 


|(2A (1) B (1) ) 


(2.131) 


with 


A - 


.4 1 «i (r o> . 9 

<o a *3 »l <r o ) + 


3r 


9 ®o (r~ ) 
IT ®3 (r 3 ) llr„ } 


(2.132a) 


B 


( 1 ) 


' 1 " 


2»j(r 0 ) 


2 ®3 <r 3 ) ,8 


<f «l (r o> + «3 (r 3> + -?T^- + f °> • 


_ -I t f(r ) a 

- <3 0i( r „) + *3 “ *?». + 11 r 3®3 r 3 


r r 


(2.132b) 

1 


( 1 ) 


3 w l'*o' 3 w l^o' 3r ( 

9 ,// 9 ®3 (r 3 ) 8 

“ *3 + lir, “ 3 CT 


2 x 2 W 


(2.132c) 


O j w / \ . £ . // / \ « 

X W + 3 "3 

2 *3 (p 3* , 16 


r +r 3®3 (r 3 } 


+ 2®3 (i*3 ) r 


+ x a 


(2. 13 2d) 
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At this point it is appropriate to make a few com- 
ments. The origin of the negative mode gammas can be under- 
stood from Eqs. (2.116) and (2.119). Note that the nearest 
neighbor third derivative of the potential is not present 
in Eq. (2.116) but is ir (2.119). Since ®i'(r o ) is 

of opposite sign to «^(r Q ) for any potential that gets 
'’softer" with increasing distance, the cancellation of 
*J(r ) in the TAX mode makes this mode gamma extremely 
likely to be negative. The cancellation of 0^'(r o ) is a 
direct consequence of the fact that the diamond structure 
has a nearest neighbor central force instability and the 
TAX mode is unstable in this case. As we shall see, the 
modes where cancellation of ®^(r o ) occurs are the modes 
that exhibit negative mode gammas. Although the C 44 
"mode" also has a nearest neighbor central force instabi- 
lity, complete cancellation of $^(r o ) does not occur in 
the mode gamma Eq. (2.131) as it does in the unstable modes 
whose eigenvectors are determined by symmetry. 

The important point of the preceeding discussion is 
to note that crystal structures possessing nearest neighbor 
central force instabilities should be extremely likely to 
have negative mode gammas and thus are good candidates for 
the occurence of negative thermal expansion. Whether or 
not negative thermal expansion actually occurs will depend 
on the harmonic phonon spectrum as well as on how many modes 
exhibit negative mode gammas. 
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From these considerations we note that several struc- 
tures can have nearest neighbor central force instabili- 
ties. In three dimensions, the simple cubic lattice, the 
body centered cubic lattice, the hexagonal diamond struc- 
ture and the graphite structure all have this feature. In 
two dimensions, the square lattice and triangular coordin- 
ation or two dimensional "graphite" both exhibit nearest 
neighbor central force Instability. It is also possible 
to obtain an instability in one dimensional chains of atoms 
that are allowed to have transverse degrees of freedom. 


rilOl DIRECTION 

As noted in Chapter 1, the mode is the "slow" 

TA mode in the elastic region. Using Eqs. (1.213) and 
(1.214) we obtain in the elastic region that 

1) 1 

3a(C ll” C 12^ 

(2.133) 



♦ 2i,'< 1 W' 1 W (1 W 1) J , 


or in terms of potential derivatives 

1 


v leD 

^4 


3a ^ll~Ci2 ) 


®,'(r ) 

•l (r o > ? + r 2®2 (r 2 ) 


, » X 7 ®2 (r 2 ) , 32 
+ 70 2 (r 2 ) — 


^ fit f v 

+ H r 3®3 (P 3* 


89 . , , , 89 *3 (r 3 ) , x 

+ u ®3 (r 3 -it ”i^r“ +4a ®4 (a) 


4® ' (a) 


+ 4$ ^ (a) 1 6c 


(2.134) 


Note the cancellation of ®J['(r o ) in the above expression. 
This does not occur for any other modes along symmetry 
directions in the elastic region. It is also interesting 
to note that experimentally, the slow TA[110] mode is 
the only mode along a high symmetry direction that exhibits 

a negative mode gamma in the elastic region for the diamond 
19 

structure . 
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In the course of this work it was necessary to have an 
analytic expression for the mode gamma of the mode at 

the K point. Although the expression Is somewhat lengthy. 
It was found very useful In obtaining a reasonable fit to 
the available experimental data. Thus we have 



* — [A (1) -(a 2 +8C 2 )"^(aa (1) +8QD (1) ) ] 

12ur 


(2.135) 


where 


J* - [?a -fc^+sr 2 )*] , 

l:k 


(2.136) 


and the definitions 


M t 


8 a- ( 2-K/2 ) g + (^ + 72 ) a+ ( 1 4+6^2 ) ix 
+ 2i'-+.(6+2«/2 ) X+24u ' + ( 4S+5//2 )y * 


a 2 

4 tp- 

+ 6^"+2X" + — (-14.706683) 


4irp l p 2 

1 A 

£ w 

s o 


r^r 

s o 


. 2 
4ttp 2 


(40.512545) + (-18.532303) 


tjt 

s o 


J (2.137) 
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M A 


u> . 


o 

_ ®9 (r 2 ) 

+ (843</2)r 2 * 2 (r 2 )+(12+5V2)(#J(r 2 )--~“ ) 

2 

+ usfcfli r3#3 - <r3 , + ii^ (t3(r3) 

60^ (a) _ 

+ 2a®^(a)+6$^(a) J (±^ + 2^5)? 


- 7 ^-[(14.706683)p 1 p{ 1) + (18.532303)p 2 p^ 1) ] 
'SO 


+ T^-( 40. 512545) (p 1 p^ 1) +p 2 p{ 1) ) 
■ s o 


(2.138) 


M A 


* 

(4-2v / 2)a-(2W2)g + (12 -i^)a-(2+2V5)n 


+ 2^ + (2+Z/2)X + (4-6«/2)^'+(14-llv / 2)l'' 

A 2 

4irp, 

- 2n"+2X"+ — ^(-17 . 790778 ) 


s o 


4irPiP 2 


4 


tr- — -(26 . 584681 ) + r7 r(- 12.350976) 


s o 


e rr 

s o 


(2.139) 
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M a (1) - P(2/3-^)r o 0^(r o ) + l 10 " 3 ^) (0 ' (r Q ) 

$X(r«) 

- 75 r 2 ® 2 fr 2 )+v ^ ®2^ r 2*“ v ^ — r 

2 

+ hm. r3#s . <rj) + afigaffli(,* (r3 , -!a£i 


2® T ( a ) *»a/5 

+ 2a®"(a)-2®"(a) + — | + (^4p “ 


24 ) a 


8 -J- t (17.790778)p 1 p^ 1) + (12.350976)p 2 p^ 1) ] 


s o 

7%-( 26 . 584681 ) ( PjPg* ^Pi* * > 
so 


(2.140) 


M C 


- 75 e a 4. (VS-8)!/' 

4irpJ 4irp..p 9 

+ ^-^(-5. 6453821) + ^-^(8.5753239) 

£ •* 6 W 

8 0 SO 

. 2 

4irp„ 

+ — ^(-5.1730388) 

€ “ 
so 


M C 


(1) 


75 


®,'(r ) 


(2.140) 


- *f(r ®:'(r o )-®:(r ) +-i— i 2 -) 

3 oio i o r Q 

. (75-8 ) /_ v ,®3 (r 3 ) v 4(1+372) 

+ ^ Ij — (r.®_ (r 0 )-® 0 (r 0 ) + — ) k 


11 V *3 W 3 V *3' *3 W 3 


- ~g— T (5.6453821)p 1 p{ 1) +(5.1730388)p 2 p^ 1) ] 
"s o A 

4 -g-(8.575323 9)(p 1 P^ 1) +p 2 p{ 1) ) 


s o 


(2.142) 
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Note that there la no obvious cancellation of «T(r ) in 

1 o 

Eqs. (2.135-2.142). However for silicon we do find a 
strong reduction of the nearest neighbor third derivative 
contribution to the mode gamm* relative to the nearest 
neighbor first and second derivative contributions. 


fllll DIRECTION 

At the L point we find from Eq. (1.237a) 


LAL 


3m “Ul 




+ 2*- (1) -46-< l W (1) 2X- (1) + f a (1) 

+ 7 ^(3p 1 -2p 2 )<3p' 1) -2p^ 1) )(tJ*(L) + 

■ o 


+ vJJ ( L>) 


or in terms of potential derivatives 


(2.143) 


LAL 




**o * 8 . 8 W 

T ®l (r o ) + 3 *l (r o ) 3 ;■ 

o 

2®;(r 9 ) 

+ 4r 2® 2 ( r 2 )’ k 2® g (r 2 ) — 

r 3 - , 32 . # , ^ 32 ®3 (r 3 ) 

+ TT ®3 (r 3 } + TT ®3 (r 3 ) "IT r— 

- .»>, x . #, V . ®4 (a) 64 

+ 2a® 4 (a)+4«, 4 (a)-4 — j- a 


+ 7 ^-(3p 1 -2p 2 )(3p{ 1) -2p^ 1) )(T^(L) + 
so 


•t vJJ(L)) 


(2.144) 


Similarly from Eq. (1.237b) 
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Y LOL 

or 

Y LOL 


For the 
Y TAL " 


1 

2 — 

3Uu ioL 


r a (1) +2e (1) +4u (1) -»-4V (1) +2X (1) ^' (1) +3X' (1 ^ 

-2 V ' (1) + 46' (1) + 4^ (1) ^X' (1 > 

+ 7T5-(p 1 +2p 2 )(Pl 1) + 2 P2 1>) < T ?? (L >- v ii a > ) 


’■ o 


(2.145) 


1 

3M “Ll 


Vl (r o )+4r 2®2 (r 2 )+2 ®2 (r 2 ) ” 


2®2 (r 2 } 


56 ®3 (r 3 } 


43 , «< / \ ,56 . » / _ \ %j 

+ TT r 3*3 (r 3 ) + TT W -IT — 

4$'(a) 

2a©"(a) + 4® * (a) j 


+ (p^+2pg )(Pj^^+2p2 ^)(T^^(L) -V^j(D) 

(2.146) 


U s o 


TA mode we obtain using Eq. (1.241a) that 


1 

3M “tal 


s < 1 >-B (1 W 1) -2p (1) + 2X (1 W (1) 

+ 3X' <1, +l/' <1) -26 ' (1) +4u' (1) +2\' <1) +3o <1> 
+ 7%-(Pl-P2 )( Pl 1) -P2 1>) <V i3 (L) - Ti 12 (L,) 


s o 


(2.147) 


or in terras of potential derivatives 


180 


TAL 


1 

3 “4al 


»i (r o ) •“^- £> - + r 2®2 <r 2 >+5 ®2 <r 2 ) 
o 

. 0 2 (r 2 ) .28 ^ , 71 , . 

-5 __ + IT r 3 * 3 (r 3 ) + TT * 3 (r 3 ) 

(r.j ) 

71 3 -j“— + 2a*"(a) + 4*'(a) 


n 


it 


-4 -JL — - 6a 


nH Pi-p ? > <P i 1) -P2 l> > <vJJ<l)-tJJ<l) ) 

s o 


(2.148) 


For the TO mode obtain using Eq. (1.241b) that 


TOL 


1 

3m 4.l 


3 . a V 1 W 1 , -*< 1 W 1 W (1> 


+ X' <li -P' ( 1 , + 26 - ( 1 W' (I) + 2 X' < 1 ) + n 0 (l)j 


2 ir 


J-(3 Pl +p 2 ) (3pj 1) +p^ 1) ) (Tj 3 (L)+V^j(L) ) 

■ 

(2.149) 


s o 


or in terms of potential derivatives 



y ?ol 


3Vu3 tol 


I *l (r o ) "I ' V ~ ' 

o 

„ , ®2 (r 2 ) 

+ r 2 » 2 (r 2 ) + 5,'(r 2 )-S -2jJ- 

, 16 _ , . 17 . 17 ®3 <r 3 > 

+ IT Vs 1 ^ 1 + TT ®3 (r 3 ) "IT ' r 3 • 

+ 2a®"(»)+4«'(») -4 a' 1 ’ 


2ir , 


s o 


( 3 Pl +P 2 ) (3 p < 1 ) +p £ 1 ) ) (Tjj (L ) +V^2 (L ) ) 

(2.150) 


At this point it is interesting to note that $'"(r ) 

l o 

is absent from the TAL mode g!,mm&. Experimentally, this 
mode also has a negative mode gamma. 


ANHARMON IC FIT TO EXPERIMENTAL DATA FOR SILICON 


In this section the method used to determine the an- 
harmonic parameters for silicon is described. Using these 
parameters we then present dispersion curves for the 
Gr^neisen parameters along symmetry directions. 

Before proceeding to the enharmonic model used in the 

present work, we want to discuss two previous enharmonic 

models used for silicon. The pioneering work of Dolling 
20 

and Cowley modeled the anharmonic potential energy by a 
two-body central potential interaction between nearest- 
neighbors. It was then assumed that two independent para- 
meters could describe the anharmonic interaction and these 
two parameters were adjusted to give a reasonable fit to 
the experimentally measured thermal expansion. However, 
the expression used to compute the mode gammas, Eq. (4.5) 
of their paper, is incorrect. Their model is also incon- 
sistent in using a shell model for the harmonic properties 
and a rigid ion model for the anharmonic properties with 
no connection between the two. Furthermore, it can be 
shown that nearest-neighbor anharmonicity of any type is in 
major disagreement with the experimentally measured mode 
gammas . 

21 

The work of Jex was along similar lines to that of 
Dolling and Cowley but Jex utilized first ind second neigh- 
bor central potential anharmonicity. Unfortunately, Jex 
made the approximation of neglecting the first and second 
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derivatives of the potential. This approximation has 
severe consequences for the mode gammas as can be seen 
from Eqs. (2.116), (2.134) and (2.148). One of the con- 
sequences is that 18 required to be positive in 

order to obtain negative mode gammas. Another consequence 
is that several mode gammas disagree with experiment by 
more than a factor of two. Jex concluded that ® 2 ^ r 2 ^ 
was more than a thousand times smaller in magnitude than 
<ft"'(r Q ) , thus seemingly justifying his assumption of only 
first and second neighbor enharmonic ity . However, using 
our analytic expressions with the approximation of Jex, it 
is easy to show that, the values of the mode gammas obtain- 
ed by Jex are inconsistent with his small value of # 2 ^ 2 ^ • 
Based on the above considerations, we seriously question 
the l.'.newidth calculation of Jex as well. Furthermore, 
one can show that more general second neighbor anharmoni- 
city is unable to give reasonable agreement with the ex- 
perimentally measured mode gammas. 

In the present calculation, we have used the ten para- 
meter harmonic model of Chapter 1 together with six inde- 
pendent parameters for the enharmonic model. The six in- 
dependent parameters are 0^(r o ) , $2 (r V ’ ®3^ r 3^ » * 

and • These parameters were obtained by per- 

forming a weighted, linear, least square fit to the fol- 

lowing experimental data Yra , y tax , y tox , Y TAL , Y T0L , 

(el) (el) (el) 

Y LA[100] • Y TA[100] * Y Eg” K an(i Y £ 4 • The experimental 


values of these mode gammas are listed in Table 5. In 
Table 6 we list the values of the enharmonic parameters 
determined from this fitting procedure. Note the rapid 
fall off of the third order potential derivatives by the 
fourth neighbor distance, though not as rapid as the fall 
off of the second order potential derivatives in Table 4. 
In Fig. 10 we have plotted the mode gammas along symmetry 
directions for this model. These curves are a marked im- 
provement over those of Jex. Though the agreement with 
experiment is at worst 35^, most of the modes are much 
better and considering the large uncertainties in the 
experimental data this is not too bad. Note that the mode 
gammas reflect the same difficulty found in the harmonic 
dispersion curves, namely, the uxpole model produces con- 
siderably more lowering of the mode gamma at the X point 
than at the L point. This causes the mode gamma at L 
to be less negative than experiment and the mode gamma at 
X to be more negative than experiment. 
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TABLE 5 

Experimental values of mode Grtineisen parameters for 
22 

silicon . Elastic mode gammas determined from the pressure 
dependence of the elastic constants. 


Yj^ " .98 ± .06 

y TAX " " 1 * 4 * * 3 
Y T0X “ 1,5 * - 1 

y tal " " 1 * 3 * ,3 
y tol " 1,3 * * 2 


^ 4 K 


1.0 ± .3 



K 


Y 


(el) 

LA[100] 


.3 ± .1 


1.11 


Y 


(el) 

TA[ 100 ] 


.325 



-.1 
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TABLE 6 

Anharmonic parameters for silicon determined by the 
method described in the text. p^ 1 ^- z i^ e and P<>~ ^ 
where e is the magnitude of the electron charge in C.G.S. 
units. 


$"'(r 0 ) - - 4.5272xlQ 13 dyne/cm 2 

0 2 ( r 2 ) " " 2 * 3833xl ° 12 
( r 3 > " 6.8387X10 11 

0 ^(a) - - 4.94to8xl0 11 

Z^ X) - - 0.405426 z^ 1) - 0.608139 
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Figure 10. Mode gammas for silicon along symmetry 
directions using the model described in the text (solid 
lines). Circles are experimental data. 









CALCULATION OF THERMAL EXPANSION FOR SILICON 


Having obtained the model parameters, we can use Eq. 
(2.59) to calculate the thermal expansion coefficient. 

The sum over the first Brlllouin zone of the crystal is 
the major difficulty in carrying out this procedure. Sev- 
eral numerical methods have been developed to evaluate such 
23 

sums. In the interest of maximum computing efficiency, 

we have chosen the special k point method of Chad! and 
24 

Cohen. Briefly, this method evaluates sums over the 
Brlllouin zone by 

BZ n 

5 f UO - N L a . f (k . ) , (2.151) 

k i-1 1 1 


with 


n 

S a. - 1 . (2.152) 

i-1 1 

Here f(k) is a smoothly varying function of wavevector, 

N is the number of unit cells of the crystal, k i are 
the special it points, the associated weighting fac- 

tors and n is the number of special it points used. As 
the number of special it points increases, the approxima- 
tion, Eq. (2.151), Improves. In the present study we have 
investigated the case of 10 special k points and 60 spec- 
ial k points. These special points are all in the irreduci- 
ble l/48th sector of the Brlllouin zone and are equivalent 
to 256 points and 2048 points respectively in the full zone. 
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The 60 special It points and associated weights for the 
fee lattice were obtained from D.J. Chadi 25 and are listed 
in Table 7 for reference purposes. 

A computer program was written to compute the fre- 
quencies and mode gammas for a general it point. The ther- 
mal expansion was calculated as well as the thermal strain 
€ given by 

BZ 

€ ■ op% - -£ w(icj)Y (itj)(n +i) . (2.153) 

3B o V o kj ltj 2 

The results of these calculations for 10 special it points 
are presented in Table 8 and for 60 special it points in 
Table 9. The thermal expansion coefficient is plotted in 
Figs. 11 and 12 for the calculation using 60 it points. 

As can be seen from Figs. 11 and 12, the agreement 
with experiment is good from 17K to 700K . The devia- 
tions below 17K are attributed to two sources. First, 
the experimental points at low temperatures are subject to 
large uncertainties due to the fact that the thermal ex- 
pansion coefficient is getting extremely small. Second, we 
suspect that the inclusion of more than 60 special k 
points will raise the theoretical curve below 17K . We 
suspect this because the calculation with 10 k points was 
raised considerably below 17K by going to 60 k points. 
The deviation from experiment near 70K is attributed to 
the fact that our fit to the TAL mode gamma is not nega- 
tive enough. As can be seen, theory and experiment are in 


1 


191 

good agreement to 700K , at vhlch point the two depart. 

28 1 

Note that the Debye temperature for ell Icon la 645K and \ 

29 

the melting point 1687K. However, even at 1600K , ] 

"a 

the difference between theory and experiment la only 12^. J 

Considering the approximate fit of our mode gammas we con- ! 

aider this good agreement. 

j 

In our calcula' -ons to this point, the bulk modulus 

i 

was taken as temperature Independent. A softening of the 
bulk modulus with increasing temperature would raise the 
theoretical curve and improve agreement with experiment. 

We can estimate this effect by using the quasihai monic bulk 

j 

modulus given by 

B - B + cB (1) , (2.154) 

with 

40 (1) -cr (1) +16l/ fl) -8X (1) -llu ' (1) 

4- 132v'^^ + 8X" (1) - 16n" (1) 

(2.155) 

and B given by Eq. (1.217). The dashed curve in Fig. 12 
is obtained by multiplying the solid curve by the correc- 
tion factor B(300K)/B(T) . Note the improved agreement 
that results from taking this temperature dependence of 
the bulk modulus into account so that now the theory and 
experiment are within 5^ at 1600K . However, it is not 


U) - T~ 
3a 


I 
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unreasonable to suggest that at these high temperatures 
higher order enharmonic effects 'could be playing a role or 
even other contributions to the free energy, such as the 
thermal creation of imperfections. It appears that these 
effects are quite small for silicon. 

As a final remark, we can combine Eqs. (2.53) and 
(2.54) to obtain 


P„ - 3B € , (2.156) 

o o o 

where B Q is the bare bulk modulus, e Q the zero- point 
strain and P Q the zero-point phonon pressure. Using e Q 
from Table 9 and from Table 3 

B q - 9.5847xl0 11 dyne/cm 2 , (2.157) 

we obtain 

P - 4.433 kbar . (2.158) 

o 

As mentioned earlier, P Q is the pressure required 
to make the strain zero at T - 0 . It is a purely quan- 
tum mechanical effect and in a result of allowing the atoms 
to move in an anharmonic potential. ThiB should be taken 
into account in calculations of the phase transition pres- 
sure of silicon from the diamond structure to the 0-tin 
structure. Such calculations are performed with the atoms 
motionless. Although zero- point corrections to the total 
energy are made, to the author's knowledge, this is the 
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first calculation of ths aizs of ths zero-point contribu- 
tion to ths prassurs. Not# that ths calculations of Yin 
and Cohen 30 gave a phass transition prsssurs of 99kbar 
and ths sxpsriasntal valus is 125kbar. Ths zsro- point 
pressure of 4.4kbar is interesting in light of this 
difference. It is an additional pressure that '•s required 
to be applied in order to have the same volume for a given 
pressure as if the atoms were all at rest, thus tending to 
raise the calculated transition pressure. However, to be 
conclusive, one should also calculate the zero-point phonon 
pressure for silicon in the 3-tin structure as well. 
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TABLE 7 

60 special lc points 
ol ( 3 ^)““ • (Courtesy D.J. 

( 1, 1,1) , ( 3, 1,1) , ( 5 
(11, 1,1) , (13, 1,1) , (15 
( 7, 3,1) , ( 9, 3,1) , (11 

( 5, 5,1) , ( 7, 5,1) , ( 9 

(15, 5,1) , ( 7, 7,1) , ( 9 

7,1) , ( 9, 9,1) , (11 

( 3, 3,3) , ( 5, 3,3) , ( 7, 

(13, 3,3) , (15, 3,3) , ( 5, 

(11, 5,3) , (13, 5,3) , (15, 

(11, 7,3) , <13, 7,3) , ( 9, 

( 7, 5,5) , ( 9, 5,5) , (11, 

( 9, 7,5) , (11, 7,5) , ( 7, 

Associated weights 

(x,x,x) points 
(x,y,y) points 
(x,y,z) points 


or the fee lattice in units 
Chadi) 

1,1) , ( 7, 1,1) , ( 9, 1,1) , 

1.1) , ( 3, 3,1) , ( 5, 3,1) , 

3.1) , (13, 3,1) , (15, 3,1) , 

5.1) , (11, 5,1) , (13, 5,1) , 

7.1) , (11, 7,1) , (13, 7,1) , 

9.1) , (13, 9,1) , (11,11,1) , 

3.3) , ( 9, 3,3) , (11, 3,3) , 

5.3) , ( 7, 5,3) , ( 9, 5,3) , 

5.3) , ( 7, 7,3) , ( 9, 7,3) , 

9.3) , (11, 9,3) , ( 5, 5,5) , 

5,5) , (13, 5,5) , ( 7, 7,5) , 

7,7) , ( 9, 7,7) , ( 9, 9,5) , 

a - 1/256 
a - 3/256 
a m 6/256 
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TABLE 8 

Volume thermal expansion coefficient and thermal strain 
for silicon using 10 special 2 points 


T 

m 

“-1 
(K X ) 

€ 3 

Cxio 3 ) 

0 

0 

1.54173 

5 

-1.81 73 4x1 O’ 13 

H 

6 

-1 . 94981xl0~ 12 

ft 

7 

-1 . 013 OlxlO” 11 

ft 

8 

-3.37142X10’ 11 

If 

9 

-8.40683X10’ 11 

f f 

10 

-1 . 73045xl0~ 10 

ff 

11 

-3 . 13907xl0~ 10 


12 

-5 . 25860xl0~ 10 

ff 

13 

-8 . 40982xl0~ 10 

M 

14 

-1.31147X10 -9 

ff 

15 

-2.01518xl0 -9 

If 

16 

-3 . 05801xl0 -9 

f » 

17 

-4 . 57273xl0~ 9 

ff 

18 

-6.71433xl0’ 9 

1.54172 

19 

-9.6528 xlO” 9 

ff 

20 

-1.35649x10’® 

f f 

21 

-1. 86243x10"® 

1.54171 

22 

-2.49942x10’® 

1.54171 

23 

-3.28196x10’® 

1.54170 

24 

-4.22219x10 

1.54168 
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TABLE 8 (cont . ) 



-1 

. (K X ) 

3 

(Wl 

25 

-5.32952x10"® 

1.54167 

30 

-1.34971xl0" 7 

1.54152 

35 

-2.56732xl0" 7 

1.54120 

40 

-4. 053l0xl0" 7 

1.54065 

45 

-5.64063xl0" 7 

1.53984 

50 

-7 . 17923xl0" 7 

1.53877 

55 

-8 . 55092xl0" 7 

1.53 745 

60 

-9.67052xl0" 7 

1.53593 

65 

-1.048 04x1 0" 6 

1.53425 

70 

- 1.09 45 Oxl 0" 6 

1.53246 

75 

-1.10464kl0" 6 

1.53062 

80 

-1 . 07810x10"® 

1.52880 

85 

-1 .01569xl0" 6 

1.52705 

90 

-9.19105x10"'* 

1.52543 

95 

-7.90721\10“ 7 

1.52400 

100 

-6 . 33396xl0" 7 

1.52281 

110 

-2. 44648xl0 -7 

1.52132 

120 

2 . 21055xl0 -7 

1.52126 

13 0 

7 . 38850xl0" 7 

1.52285 

140 

1 . 28710xl0" 6 

1.52622 

150 

1 . 84817xl0~ 6 

1.53145 

160 

2 . 40845xl0~ 6 

1.53854 

170 

2.95791xl0" 6 

1.54749 
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TABLE 8 (cont . ) 



T 

W 

a _l 

(K A ) 

€ 3 
11 

180 

3.48954xl0“ 6 

1.55824 

190 

3.99872xl0“ 6 

1.57073 

200 

4. 48 268x1 O*" 6 

1.58487 

210 

4.94002xl0~ 6 

1.60059 

220 

5. 3 7036x1 0~ 6 

1.61778 

230 

5 . 77401xl0~ 6 

1.63636 

240 

6 . 15176xl0 -6 

1.65624 

250 

6.50470xl0~ 6 

1.67734 

260 

6.83412xl0~ 6 

1.69958 

270 

7. 1413 9x1 0" 6 

1.72288 

280 

7. 42792x1 0“ 6 

1.74717 

. 290 

7 . 69511xl0~ 6 

1.77238 

300 

7.94432xl0" 6 

1.79845 

350 

8 . 96355xl0 -6 

1.93983 

400 

9 . 69659xl0~ 6 

2.09565 

450 

1 . 02353xl0** 5 

2.26197 

500 

1.06401xl0” 5 

2.43609 

550 

1.09505xl0“ 5 

2.61612 

GOO 

1 . 1193 2xl0~ 5 

2.80073 

650 

1 . 13860xl0** 5 

2.98895 

700 

1 . 15415xl0” 5 

3.18006 

750 

1 . 16687xl0 -5 

3 .37351 

800 

1.17739xl0“ 5 

3.56889 
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TABLE 8 (cont . ) 



T 

(K) 

a -l 
(K X ) 

e 3 

(xlO 3 ). 

850 

1.18619xl0 -5 

3.76588 

900 

1 . 19361x10”® 

3.96421 

950 

1.19994xl0~ 5 

4.16369 

1000 

1.20537xl0~ 5 

4.36414 

1050 

1.21006xl0“ 5 

4.56544 

1100 

1.21414xl0“ 5 

4.76746 

1150 

1.21772xl0~ 5 

4.97012 

1200 

1.22087xl0“ 5 

5.17334 

1250 

1 . 22365xl0” 5 

5.37706 

1300 

1.22613xl0 -5 

5.58121 

1350 

1 . 2283 4xl0 -5 

5.78575 

1400 

1 . 23 032xl0” 5 

5.99064 

1450 

1.23210xl0 -5 

6.19585 

1500 

1.23371xl0" 5 

6.40133 

1550 

1 . 23517xl0” 5 

6.60708 

1600 

1 . 23650xl0*" 5 

6.81305 

1650 

1 ,23770xl0" 5 

7.01924 

1700 

1 .23881X10 -5 

7.22561 
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TABLE 9 

Volume thermal expansion coefficient and thermal strain 
for silicon using 60 special lc points. 

T a -l € 3 

00 (* ) (xio 3 ) 


0 0 1.54184 


5 

1 . 41 794x1 0" 1 ^' 

ff 

6 

4.28023X10” 11 

ff 

7 

9.S5599X10 -11 

M 

8 

1.86983xl0" 10 

ft 

9 

3 . 28359xl0** 10 

ff 

10 

5.31819xl0" 10 

ff 

11 

7.97764xl0” 10 

ff 

12 

1 . 10894xl0” 9 

ff 

13 

1 . 424472xl0” 9 

f f 

14 

1.67673xl0“ 9 

f f 

15 

1.76651xl0“ 9 

ff 

16 

1 . 56576xl0 -9 

ff 

17 

9 . 19612x10”*° 

ff 

18 

-3 . 47348xl0~ 10 

ff 

19 

-2. 423 03x1 O’ 9 

f f 

20 

-5 . 49874x10" 9 

ff 

21 

-9. 76004x1 O -9 

ff 

22 

-1.53 785x1 0~ 8 

ff 

23 

-2.25047xl0“ 8 

1.54183 

24 

-3 . 12632xl0** 8 

1.54182 


TABLE 9 (cont.) 



T 

(K) 

H 

1 

a 5* 

€ 3 
(xlO J ) 

25 

-4.17490x10"® 

1.54181 

30 

-1 . 21275xl0~ 7 

1.54168 

35 

-2.41884xl0" 7 

1.54138 

40 

-3.89885xl0" 7 

1.54086 

45 

-5 . 48352xl0" 7 

1.54008 

50 

-7 . 02070xl0" 7 

1.53903 

55 

-8 ,39167xl0” 7 

1 .53775 

60 

-9 . 51091xl0~ 7 

1.53625 

65 

-1.03 206x1 0“ 6 

1.53459 

70 

-1 . 07851x10”® 

1.53283 

75 

-1 . 08864xl0" 6 

1.53102 

80 

-1 . 06211xl0~ 6 

1.52922 

85 

-9 . 99695xl0" 7 

1.52750 

90 

-9.03106xl0“ 7 

1.52591 

95 

-7.74721xl0" 7 

1.52450 

100 

-6 . 17395xl0~ 7 

1.52334 

110 

-2 . 28646xl0~ 7 

1.52191 

115 

-3 . 79891x10"® 

1.52171 

120 

2.37060X10” 7 

1.52190 

130 

7 . 54856xl0 -7 

1.52354 

140 

1 .30311x10"® 

1.52697 

150 

1 . 86418x10”® 

1.53224 

160 

2 . 42446xl0~ 6 

1 .53939 
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TABL)S 9 (c<»nt . ) 


(K) 

-1 

— (K X > 

* 3 

(xlO* 5 ) 

170 

2.97393x10”® 

1.54840 

180 

3.50556x10”® 

1.55920 

190 

4.01473x10"® 

1.57174 

200 

4.49869X10" 6 

1.58594 

210 

4.95604x10”® 

1 .60170 

220 

5.38638x10”® 

1.61895 

230 

5 . 79003xl0” 6 

1.63758 

240 

6 . 16777xlo” 6 

1.65752 

250 

6 . 52071x10”® 

1.67867 

260 

6 . 85013x10”® 

1 . 70096 

270 

7.l5740xi0“® 

1. 72432 

280 

7.44393x10”® 

1.74866 

290 

7.71112x10”® 

1.77392 

3 00 

7.96033x10“® 

1.80004 

350 

8.97955x10”® 

1.94169 

400 

9.71260x10“® 

2.09778 

450 

1 . 02513xl0“ 5 

2.26437 

500 

1 . 06561xl0” 5 

2.43875 

550 

1 . 09665xl0“ 5 

2.61905 

600 

1 . 12092xl0” 5 

2.80393 

650 

1 . 14020xl0” 5 

2.99241 

700 

1.15575x10”® 

3.18379 
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TABLE 9 (cont. ) 


CK) 

-1 

(K 1 ) 

(xlO 3 ) 

750 

1.16847xl0" 5 

3.37751 

800 

1.17899xl0~ 5 

3.57316 

850 

1.18779xl0" 5 

3.77041 

900 

1.1 95 21x10" 5 

3 . 96901 

950 

1 . 20154x10"® 

4.16875 

1000 

1 . 20697xl0" 5 

4.36947 

1050 

1 . 21166xl0" 5 

4.57103 

1100 

1 . 21574xl0" 5 

4.77333 

1150 

1 . 2193 2xl0 -5 

4.97625 

1200 

1 . 22247xl0" 5 

5.17974 

1250 

1 . 22525xl0" 5 

5.38372 

1300 

1 . 22773xl0 -5 

5.58814 

1350 

1 . 22994xl0" 5 

5.79295 

1400 

1 . 23192xl0" 5 

5.99811 

1450 

1 . 23370xl0 -5 

6.20358 

1500 

1 . 23531xl0" 5 

6.40933 

1550 

1 . 23677xl0" 5 

6.61534 

1600 

1 . 23809xl0" 5 

6 . 8215fi 

1650 

1.23930xl0" 5 

7.02803 

1700 

1.24041xl0‘ 5 

7.23468 
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Figure 11. Volume thermal expansion coefficient for 
silicon calculated using 60 special ic points (solid 
line). Circles are experimental data (Ref. 26). 
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Figure 12. Volume thermal expansion coefficient for 
silicon calculated using 60 special & points (solid 
line). Circles are experimental data (Ref. 27). The dash- 
ed line is calculated by applying a quasiharmonic correc- 
tion to the bulk modulus. The Debye temperature for sili- 
con is 645K and the melting point is 1687K . 
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Chapter 3 


Anomalous Thermoelastic Effect in Silicon 


"For God so loved the world, 
that he gave his only 
begotten Son, that whosoever 
believeth in him should 
not perish, but have ever- 
lasting life." 

John 3 : 16 


In this chapter the equations of visco and thermo- 
elasticity are discussed and the dispersion relation for a 
cubic material presented. The thermoelastic contribution 
to the acoustic attenuation in silicon is computed from 
1-300K . Strong attenuation anomalies associated with 
negative thermal expansion are found in the vicinity of 
17K and 125K . Comparison with experimental results is 
discussed. It is suggested that anharmonic effects are re- 
sponsible for the anomalies seen in low frequency measure- 
ments and that several materials should exhibit similar 
behaviour . 

Various studies of acoustic attenuation in silicon 
have been carried out in the past several years. Measure- 
ments of acoustic attenuation at relatively high frequencies 
have been made by Mason and Bateman^, who studied silicon 
and germanium at 300 -500 MHz and found that the attenua- 
tion went smoothly to very low values at low temperatures. 

A quantity which affects the acoustic attenuation is the 
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210 


viscosity tensor, whose elements have been determined by 
2 

Helme and Kinis for silicon at 1 GHz and 300K . On the 
other hand, measurements at low frequencies indicate 

3 

attenuation anomalies. McGuigan et_ al, measured the 
mechanical Q of a large single crystal of silicon at 20 
KHz and observed attenuation anomalies at 13K and 115K . 

4 

Davis has also observed a dip in Q near 13K at 3.4 
KHz in a large single crystal of silicon. McGuigan et al , 
attribute the anomalies to impurities and mechanical imper- 
fections of the crystal. However, no conclusive evidence 
is presented for this mechanism, and it is possible that 
the effect is an intrinsic property. It is the purpose of 
this investigation to show that anharmonic attenuation can 
very rapidly in a narrow range of temperature and exhibits 
anamalies near 13K and 125K in silicon. 

Consider the combined thermoelastic and viscoelastic 

5 . 6 

effects in a cubic crystal. The equations of motion are 



2 

a-r a 9 u v 


+ ( c 44 + ^44 8t ) 




(3.1a) 
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9 2 u. 


3 2 u 


P -JJ# " " Bor If + (C 11 + ^11 Tt * (C 44 + ^44 It* 


d 2 U d 2 U 0 

x \'TT^~TT) + (C 12 + T1 12 ^ + C 44 +r| 44 7t } x 
dx dz 


8 2 u S*u 

* ISTlf + -573t. 


(3.1b) 


3 2 u 


z _ 


at' 


a 2 

- ■« f + ‘Cu+’tu -^f + < C 44 + l44 x 


a 2 a 2 

x i tt + " 7 ^) + (C i 2 + ^12 inr +C 44 + ^44 x 


(3.1c) 



while the heat equation is 


-rp (C “C ) 1 ^ O 

c v It + -V^w <v,u) - <7 * T • 


(3.2) 


Here p is the mass density, and -n A j are the elas- 

tic moduli and viscosity tensor elements respectively, B 
is the isothermal bulk modulus, a is the coefficient of 
volume thermal expansion, T is the temperature, C p and 
C v are the heat capacity per unit volume at constant pres- 
sure and volume respectively, tc is the thermal conducti- 
vity, and the u i are the cartesian components of the 
elastic displacement. 
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We seek plane wave solutions to Eqs. (3.1) aud (3.2) 


of the form 


u. ■ u. e 
i lo 


ik* x-iort 


(3.3) 


T - T + T, e 

O 1 


ik • x- iurt 


(3.4) 


where k Is the wave vector and to is the frequency. For 
simplicity, we specialize to a longitudinal wave propaga- 
ting in the [111] direction. Substituting Eqs. (3.3) and 
(3.4) into Eqs. (3.1) and (3.2) we find that the disper- 
sion relation is 


k ± " (1-iorr ) + 1 (b + 1-iot ) 


* (1-iorr) + i(b+ l3 


4iab 1 £ 
(1-iorr)) 


(3.5) 


where we have defined k - k and 


2 2 

B a To 

2pKvf 


(3.6) 


C 11 +2C 12 + 4C 4 4 \ ~ 


(3.7) 


T 'll +2r 'l2 +4r '44 


(\l + 

\ C 11 + 


C 12 +4C 44 


(3.8) 
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The lower sign In Eq. (3.5) corresponds to t.he elastic 
mode while the upper sign corresponds to the thermal mode. 
Note that Eq. (3.5) is valid for any purely longitudinal 
wave, provided we use the appropriate velocity v Q and 
viscous relaxation time t for the direction of propaga- 
tion. Also note that t here is not the thermal relaxa- 
tion time. 

Various cases of Eq. (3.5) are of interest. For a 
pure viscoelastic effect we have 

2 J 

k - - 3 — £ . (3.9) 

v ( 1 -iu/r ) 
o 


This can be obtained from Eq. (3.5) by putting a equal to 
zero. Letting 


k - + ikg 


we obtain for urr«l 


k- - 


V 


oj 2 t 


(3.10) 


(3.11) 

(3.12) 


Note that urr«l is satisfied in silicon at room tempera- 
ture even at 1 GHz . 

For the pure thermoelastic effect t - 0 . This dis- 
persion relation has been obtained previously , 7 and in the 



case that b»a,c we have 


2 2 2 
KB a To) 

2 - 5 r 2 (3.13) 

2 P v o C v 


and la given by £q. (3.11). Equation (3.13) Is the 

usual result,* and applies to silicon at 20 KHz from 
1 -300K . 

For combined thermo and viscoelastic effects, one may 
show by expanding the dispersion relation (3.5) that in 
the case b»a,c and cjr«l , is given by Eq. (3.11) 
and 




KB 2 a 2 To. 2 

s nr" • 

2pvC 
* o v 


(3.14) 


The usual situation is that the second term of Eq. (3.14) 
is only a few percent of the first term. 

Since the temperature dependence of the viscosity ten- 
sor elements has not been measured at low frequency over a 
range of temperature and has not been calculated using a 
realistic phonon model for silicon, we focus on the ther- 
moelastic contribution to the attenuation. Using data 

8—10 

found in the literature, we have evaluated Eq. (3.13) 
for silicon in the range of 1-300K at 20 KHz . These 
results are presented in Table 10. Figure 13 is a plot of 
these results. Figure 14 is an expanded version of the 
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region from 1-50K . Note that the plots are normalized 
to be independent of frequency in the range of validity of 
Eq. (3.13). 

The attenuation shows rather marked features. The 
maximum at 14K is associated with the thermal expansion 
coefficient passing through a maximum near 14K . The dip 
near 17K is associated with the thermal expansion coeffi- 
cient going through zero and changing sign. The broad peak 
from 35 - 60K is due to the combined factors of the 
thermal conductivity passing through its maximum and de- 
creasing, while the magnitude of the expansion coefficient 
increases. The slight shoulder in the curve from 70 - 80K 
is due to the thermal expansion coefficient passing through 
its minimum negative value. The dip in attenuation near 
125K is due to the thermal expansion coefficient changing 
sign from negative to positive in this region. Note that 
the maximum thermoelastic attenuation occurs in the vici- 
nity of 40 - 45X and not at room temperature. In fact, 
the thermoelastic attenuation at 40K is more than a 
factor of 3 times greater than it is at room temperature. 

In this region we expect the thermoelastic effect to be 
more than a few percent of the total attenuation for longi- 
tudinal waves. Careful measurements of the viscosity ten- 
sor elements at low frequency are necessary to substantiate 
this conclusion. 

The existence of thermoelastic attenuation anomalies 


3 4 

near the temperatures at which McGuigan et al. and Davis 

observed anomalies in the mechanical Q suggests that an- 

harmonic effects are the cause of the latter anomalies and 

not impurities or mechanical imperfections in the crystal. 

Indeed, the sample of McGuigan et, al. was boron doped with 

15 3 

a concentration of 3 x 10 boron atoms/cm , while the 

14 

sample of Davis had a boron concentration of 8 x 10 

3 

atoms/cm , both of which are rather low impurity concen- 
trations. In order to settle the question, it is necessary 
to perform experiments on extremely pure samples and see 
if the anomalies persist. If they do, one would have very 
strong evidence in favor of the anharmonicity mechanism. 
Furthermore, since several semiconductors such as Ge , 

GaAs , ZnS , ZnSe and CdTe exhibit negative thermal expan- 
sion,^ one may expect that they also will have anomalous 
acoustic attenuation if, in fact, the effect is intrinsic. 
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TABLE 10 

Values of kg for silicon at 20 KHz along the 
[111] direction using Eq. (3.13). Units of kg are 
10” 12 /meter . 


Y(K) 

k 2 

T(K) 

k 2 

T(K) 

k 2 

1 

.01164 

20 

47.87 

135 

17.58 

2 

.1644 

22 

181.1 

140 

31.42 

3 

.7736 

25 

510.3 

150 

62.62 

4 

2.324 

30 

1049 

160 

104.3 

5 

5.408 

35 

1365 

170 

147.2 

6 

10.38 

40 

1467 

180 

188.2 

7 

17.24 

45 

1467 

190 

223.6 

8 

26.13 

50 

1387 

200 

251.8 

9 

36.78 

60 

1141 

210 

288.4 

10 

48.48 

70 

555.7 

220 

319.6 

11 

60.17 

80 

507.9 

230 

344.0 

12 

72.17 

90 

266.4 

240 

364.0 

13 

87.42 

100 

113 .7 

250 

3 78.4 

14 

94.11 

1C5 

65.46 

260 

398.0 

15 

62.08 

110 

33.83 

270 

414.6 

16 

30.64 

115 

12.29 

280 

429.3 

17 

. 06894 

120 

1.883 

290 

438.0 

18 

.8661 

125 

.4296 

300 

448.1 

19 

16.99 

130 

6.347 
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Figure 13. Thermoelastic attenuation for silicon 
along [111] direction from Eq. (3.13). 
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Figure 14. Thermoelastic attenuation for silicon 
along [111] direction from Eq. (3.13). 


1 
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Chapter 4 


Conclus ions 


"And this is the record, 
that God hath given to us 
eternal life, and this life 
is in his Son. He that 
hath the Son hath life; and 
he that hath not the Son of 
God hath not life." 

I John 5:11,12 


In this final chapter we summarize the main features 
of the work presented in the previous chapters and offer 
some suggestions for future research. 

It has been found that long-range, nonlocal dipole- 
dipole interactions are able to provide lowering of the 
frequencies of the TA modes in diamond structure crystals. 
The dipoles also have the feature of not affecting the 
elastic constants or Raman frequency. It was necessary to 
include central potential interactions to fourth neighbors 
in order to fine tune the agreement with experiment. A 
reasonable fit to the experimental data has been obtained 
that exhibits a rapid fall off of the first and second order 
potential derivatives with increasing distance. 

It is believed that further investigation of long- 
range electrostatic interactions in diamond structure 
crystals will be profitable. Specifically, one might 
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examine the effect of including the local quadrupoles of 
1 

Lax, in addition to the nonlocal dipoles, in an effort to 
obtain agreement with experiment using a smaller number of 
parameters. The resulting dipole-quadrupole and quadrupole- 
quadrupole interactions should improve the problem of the 
angular variation of the dipolar Interaction energy. In 
addition, one might also examine the effect of the short- 
range corrections to the dipole-dipole interaction energy 
as discussed in Appendix B. These short-range corrections 
have a different angular variation than the usual dipole- 
dipole interaction energy. 

An anharmonic model has been developed that is a con- 
sistent extension of the harmonic model. Using this model 
we have obtained analytic expressions for the mode Gruneisen 
parameters that explain in a simple way the origin of the 
negative mode gammas. The relation of negative mode gammas 
to nearest neighbor central force instability is illustrat- 
ed by these analytic expressions. It has been found that 
the approximation of neglecting first and second order po- 
tential derivatives in comparison with third order potent- 
ial derivatives is particularly severe in the diamond 
structure. It has nlso been found that long-range anhar- 
monicity is necessary to explain the experimental mode 

gammas in silicon and that the short-range anharmonicity 

2 3 

calculations of Dolling and Cowley and Jex are inadequate 
to obtain even moderate agreement with experiment. Using 
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the model developed in the present work, a reasonable fit 
to the experimental mode gammas for silicon has been ob- 
tained that exhibits a rapid fall off of the third order 
potential derivatives with increasing distance. 

The volume thermal expansion coefficient of silicon 
has been calculated from 5K to 1700K using the above 
lattice dynamical model. Surprisingly good agreement with 
experiment from 17K to the melting point of 1687K has 
been obtained. This shows that higher order anharmonic 
contributions to the thermal expansion for silicon are very 
small, even at temperatures near the melting point. To fur- 
ther improve agreement with experiment below 17K , it is 
suggested that more than 60 special lc points be used in 
evaluating the sum over the Brillouin zone. It is also 
believed that a better fit to the experimental mode gamma 
at the L point would improve agreement with experiment 
in the region of 75K , where the thermal expansion be- 
comes the most negative. 

Associated with the negative thermal expansion, it 
was found that silicon has an anomalous thermoelastic 
effect. It has beer: shown that this contribution to the 
acoustic attenuation exhibits strong attenuation anomalies 
in the vicinity of 17K and 125K . This investigation 
is preliminary in nature and demonstrates that anharmonic 
attenuation can vary rapidly in a narrow range of temper- 
ature. It is the present author's conviction that a 
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calculation of the viscosity tensor for silicon as a func- 
tion of temperature will exhibit anomalies as well. How- 

4 

ever, exsistlng expressions for the viscosity tensor of 
a solid can be objected to on fundamental grounds, and be- 
fore a proper calculation can be made, a new theory of the 
viscosity tensor must be developed. The present author has 
investigated such a theory but since the results are not 
yet complete, they will not be reported on here. 
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Appendix A 


ANGLE BENDING CONTRIBUTION TO THE 
POTENTIAL ENERGY 

In this appendix we derive an expression for the 
change in angle between any three atoms in a crystal cor- 
rect to first order in the atomic displacements. The dis- 
cussion is similar to that given by Trul linger 1 and iB in- 
cluded here for the sake of completeness. 

From the geometry of Fig. 5 and the definition of the 
dot product we have 

cos e (XK |X'k '|X"k") - (A-l ) 

|x(AK|i # K'>||Sc(AK|A'g')J ’ 

where 

x(Xk|X'k') - x(lK) - x(X'k') , (A-2 ) 


and 


X(XK ) - R(XK ) + u(XK ) . 


(A-3 ) 


Now 


x(XK \Z'k')'x(Ik |X"k") 


jfctc (x'iC ') *^(iK |X" k") + 

+ $(XK |X'k')*u(XK I 1"k") 

+ £(Xk!X"k">.u(XK|X'k') + 

+ u(XK | X '< ') * u(XK J X" k *) 

m m 

(A-4) 
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where 


- ft(ZK) - $(Z'0 , (A-5) 


and 


u(Atc | Z'k ') - u(Zk) - u(Z'k') . (A-6 ) 


We also have 

|x(Zk|z'k')| 


|£(Zk|Z'k') I 


* ■ 
1 | 2^(Atc|X*K*)»u(AtCll*K*) 

I / i u> I * * * \ I 2 


| R( XK | Z'tc ') (■ 

iSCiKU'K*))* 
|it(ZK |Z'ic')p Z 


Using the expansion 


(A-7 ) 


(i+y ) _1/2 - i- |y+|y 2 +-" , 


(A-*) 


we have to first order in the displacements that 


|x(Zk|z'k') I" 1 - |R( Zk I 1 ' k ') f 1 


1 - RUk [ 111 ) 

| |C ( Z K | Z ' K ' ) | ^ . 


(A-9) 


Combining Eqs. (A-l), (A-4) and (A-9) we obtain to first 
order in the displacements that 

cos e(K|Z'K'|Z'K') - COS 6 (0) (ZK I I'k '| Z'k') - 
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- COS 0 (0) (IK | I'k' |A'k') fi(A* | A'<*)»u(AK 1 1'<') 

7TT7777*ZT7i2 


|$(Ak|A'k')| 


- COS 0 (O) (Ak|A'k'|A'k')| 


ft(AK U # K*)«u(AK I 1 " k ') 

I ^ / I I S * * V I • 


|i?(AK|A'K')! 


+ ft(AK|i*K*)«u(AKli*K*) + ftUKlA'K'VttCAKirO..-— 
|^(XK I jfc'K ') I |3 (Ak|A'k')| |3(Ak|A'k')| |i*UK|A'K')! 


(A-10) 


where we have defined the equilibrium angles by 


cos 0 (O) (Ak|a'k 'Ia'k 1 ') 


ft(jK |i*K*)»]£(iK U'k*) 

|R( Ak I A'k') I | R( Ak I A'<') I 

( A-ll ) 


Defining the change in angle from the equilibrium angle by 

A 0 (Ak|A'k'!a'k') - 0 (Ak|a'k'|A'k') - 0 (O) (Ak|a'k'|a'k') , 

(A-12 ) 


we can use the fact that 

cos 0“cos(8^^+A0)"Cos 0^°^cos A0 - sin S^^sin A0 , 

(A-13 ) 

so that to first order in A0 

cos(e (0) + A0) - cos 0 (O) -A0 Sin 9 (0) , (A-14) 


or 
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-( 


A0 - ( cot 0 


(0) cos 0 


■ In 6 


m) 


(A-15) 


Using Eq. (A-15) together with Eq. (A-10) we obtain 


A 9 (AK|A'k'| A'k') - 
- cot 0 ( °^ ( AK | A # K ' | A # K # ) 


ft(AK|A'K')»u(AK|A'K') 
; |ituK| a'k') | 2 


+ cot 9 (Q) (AK[A'K'|A*K # )| R(AKli # K # )»u(AKfA <, K # ) 


• .. • v . . 


|it( AK | A # K ' ) | 2 


- CSC 8 (Q> (AK[A'K^[A <, K J, )[ g(jftKfA^K*).u(AKfA"K") 


- CSC 8 ( 0 ) (Ak|a'k'|a'k')| r(Ak|i'k')«u(Ak|i'k ') 


|S( AK ! A'k') f |j?( AK | A*K ') 




!R(AK |A'k') I |I(AK |A'k') I 

( A- 16 ) 


Def ining 


ri(AKf A'K'fA'K') - 
- cot 8 ( 0 ) (Ak[A'k'|A'k') 

CSC e ( 0 ) (AK |a'k'|a'k') 


fl(AK | A'k' ) 
|£(Ak| a'k') | 2 

ft(AK | A' k') 


)^(AK I a'k ' ) | |R( Ak | A'k') I 


(A- 17 ) 


Eq. (A-16) can be rewritten as 

A 0 (AK | A'k ' ( A'k ') - ii(AK (A'k ' ( A'k ') » u(AK (A'k ') 

+ ti(Ak | A'k ' | A'k')»u(AK | A'k ' ) # 

( A— 18 ) 
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or in component form 

aeuic }i'<' |a'k') - £ r\ a (t< |a'k # |i'K # )u a <iic|f'K') 

a 

+ Lr\ (iKli'K'IX'tc')u (iKli'K 0 ) . 

(A-19) 

Note that Eq. (A-18) depend* explicitly on differences of 
displacements so that it satisfies infinitesimal transla- 
tion invariance. Using Eq . (A-18) it is easy to show in- 
finitesimal rotation invariance. Consider the infinitesi- 
mal rotation 

u(XK) - ux(J(i «>-3 0 ) . (A-20) 

* 

where ju>| is the infinitesimal rotation angle and 
is the origin of the rotation. Then we have 

u(!K|i'<') - wxi&CiK li'O . (A-21 ) 

Substituting Eq . (A-21) into Eq. (A-18) and using the 
vector identity that 

t< (BxC ) - -B* (Ax?) , (A-22) 


we obtain 
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aeu* \i'k' \i'k 0 ) 


-Oil 


ri(iK (1'k'|4'k')x^(1k|jI'k ') 


|x # k # |i'k')xI(Ak|A # k')J . 


(A-23 ) 


From Eq. (A-17) we see that 

SukIi'k' |I # k') X 2(1k|a'ic') - 

- esc 0 (0) (XK I Z'i; ' I x'k # ) — I (JUJ A '. K *i*ft.L* K .l *'.*') 

|S(Xk|x'k')| |$(Xk|x'k')| 

(A-24) 


and similarly 


*n(AK \l\ * (X'K ') X $(XK I jTO - 

CSC e (0> (iK |1 'k' |i'« ') . 

|$(XK | X'tC ') [ f !C( AK I i \ ') ! 

(A-25 ) 


In obtaining Eq. (A-25) we have made use of the property 
that 


e(XK ]l\ '|x'k') - Q(iK\l"<" \1 'k') . (A-26) 

Adding Eqs. (A-24) and (A-25) we obtain 

A6(Xk|X'k'|x'ic') - 0 , (A-27 ) 

under the infinitesimal rotation Eq. (A-20) . Thus the 
change in angle Eq. (A-18) satisfies infinitesimal rotation 


invariance . 
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Appendix B 


ELECTROSTATIC INTERACTION ENERGIES IN 
A NONLOCAL DIELECTRIC MEDIUM. 

In this appendix we discuss the interaction energy of 
two charge distributions in a particular type of linear, 
nonlocal dielectric medium. Results for the interaction 
energies of the first few cartesian multipoles are pre- 
sented. Further application to the case of isotropic non- 
locality is made. The special case appropriate to semi- 
conductors is examined. 

Consider the Interaction energy of two charge distri- 
butions in a linear dielectric medium. We have* 


*12 - - (B * 1) 

where 

- 4 rp 1 g * 4ir P 2 * (3-2a) 

5 - +S 2 2 - +f 2 . (B-2b ) 

Introducing the Fourier transforms 

f(x) - J*d 3 x e^^fCk) (B-3a) 

f (k) i— oJd 3 x e _i *‘ x f(x) , 'B-3b) 

(2r) J 
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we suppose for the medium of Interest that 

$(£) - €(£)$(£) . (B-3b) 

This corresponds to examining the diagonal part of the 

dielectric matrix, or neglect of umklapp terms for a cry- 
2 

stal. These terms have been shown to be small for silicon 

3 

and germanium as anticipated by Penn. After some algeb.a 
and use of Eqs. (B1-B4) we obtain 

W 12 " ^fj’d 3 x J'd 3 x'g(x-x')p 1 (x)|02^') , ('■ ■ .) 


with 


g(x) 



1®. 


iit‘x -ijk* 

“TOT 


x 



(B-5b) 


Note also that 

g(x) - g(-x) . (B-5c) 


All the nonlocal behaviour resides in the function g(x) . 

Now consider the first few terms of a cartesian multi- 
pole expansion about the point x . We have 


p(x) - Q6(x-x') -E P q 6 (x-x' ) 


a ~ a 
.2 


+ » a.t d V*P 


6 (x-x ') 


+ (higher multipoles) , 


(B-6 ) 


237 


where 6(x) 
tion of Eq. 


is the Dirac delta function. 

Direct integra- 

(B-6) gives 



Q - 

I d 3 xp(x) 
V' 

(B- 7a) 

p a ’ 

J* d 3 x x p(x) 
V' 

(B-7b) 

Q a0 " 

J* d 3 x x x B p(x) , 

v' a p 

(B-7c ) 


where V' is the volume to which the charge density is 
localized. Eq. (B-6) can be obtained by application of 
Poisson's equation to determine the charge densities 
responsible for monopole, dipole and quadrupole potentials. 
Utilizing Eq. (B-6) in Eq. (B-5a) we obtain 


W 12 


J'd 3 xj*d 3 x 'g (x— x * ) [Q x 6 (x-j^ )- 1 (1 )~-6 (x-j^ ) 


a) teVr 6 J * t V <*'-*2> 


or, 3 


a 

a p 


H M* V 

(B-8) 


Multiplying out the various terms yields 


w - w MM + w Md 4. 4- w dd 4- w d ^ 4- 

W 12 W 12 + "l2 + W 12 + W 12 + W 12 + W 12 


+ (higher raultipoles) , 


(B-9) 


238 


where 


“ Jd 3 x,fd 3 x'g(x-x')Q 1 Q 2 6(x-x 1 )&(x'-x 2 ) (B~10a) 


-Jd 3 xjd 3 x ' g ( x-x ' ) [ Q, 6 ( x-x x ) E p ( 2 )-gp- 6 ( x ' -x 2 ) 

a a 

+ Q 2 i( *'"* 2 )Sp o ( 1 ) dr’ 6 ( *“*l )J 0 


- ^J*d 3 xJ‘d 3 x < g(x-x')[Q 1 Mx-x 1 ) S Q a g( 2 ) dx ; S "-rMx'-x 2 ) 

* /w K /V H 


a, p a p 


+ q 2 6 (;'-j 2 ) o 2 


E J*d 3 xj'd 3 x'g(x-x')p (l)Pg(2) 
36 (x'-x 2 ) 1 


36 (x-x^ ) 


- - 4 E J , d 3 xJ , d 3 x'g(x-x) P^C 1 ) 

a.n.i' L 


36 (x-x^ ) 
3 x 


) 36(x'-x 2 ) 

2 M /o\ • 


3 o (x '-x 2 ) 

x v (2) ex' ex' ~ + V 2) 

p, 1 / 


36 (x-x- ) 

x V ll> ex ex 

^ pi ^ 


I E Q g(l)Q ( 2 )Jd 3 xJ , d 3 x'g(x-x') x 

ae^ 

|"/ 3 2 6 (x-x-)\ / o 2 6 (x'-x 2 )\"1 


¥)( ; 


3 x 3 x„ 
M- ^ 
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The terms in Eqs. (B-9) and (B-10) correspond to monopole- 
mo no pole interaction, monopole-dipole, monopole-quadrupole, 
dipole-dipole, dipole-quadrupole and quadrupole-quadrupole 
interactions, respectively. We now proceed to evaluate the 
various integrals in Eqs. (B-10). 

For the monopole-mono pole term Eq. (B-lOa) yields 
immediately 


W^2 * QjQ 28 »( x i" x 2 ^ • 


(B-ll ) 


Thus we see that g(x) is essentially the modified 
Coulomb's law for the nonlocal medium. 

Consider next the monopole-dipole term. We can 
immediately perform one integration, the remaining integra- 
tion to be done by parts so that Eq. (B-lOb) becomes 



_ r dgCxj-Xg) 

2 kv 2) k ~ + <*<*-«> 

a L 


dg(x 1 -x 2 ) 


2a 


‘2*a 


3x 


la 


(B-12) 


Two partial integrations give for the monopole-quadrupole 
term 


w: 


MQ 

12 


\ £ 
2 a,* 


L Q l Q aP (2) dx 


2 -> -> 2 
3 g(x 1 -x 2 ) d g(x x -x 2 ) 

2<* dx 2g Q 2 Q a3 ( ) 3x la dx ip J 


(B-13 ) 


and for the dipole-dipole term 
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w dd „ y 

W 12 u 
a,0 


»*i a »*2l> e<X l‘ X 2 ) 


P a (l)p 0 (2) . (B-14] 


Three integrations by parts give for the dipole-quadrupole 
term 


W 


dQ . 
12 


1 

2 




d g(x- -X„) 


V 1} V (2) i* 


la ax 2n dx 2i/ 


a g(x--x„) 

+ p a <2)( W 1, ex 2o dx lii dx 2l , J 


(B-15) 


Finally, the quadrupole-quadrupole terra is obtained by 
four partial integrations yielding 


W 


QQ . 1 


12 


1 L 
4 an 

v-v 


dx. 


a gUj-xg) 

dx- Q dx 


lg UA 2n dx 2i/J 


Q *e (1) V (2) 


(B-16) 


Having obtained these general results, we now examine 
the case when 


e(lc) - e(k) where k - | it J . (B-17) 

Using the definition Eq. (B-5b) and integrating over angles 
we obtain 

g(x) - , (B-18a ) 


where 
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- / v x i» sin kr 

«(r) - - J dk 


Using the form of Eq. (B-18&) and differentiating 


3 g(X--X 9 ) r . 

^i a ' Wr • [ ( w i i 2 > f( i*i-’‘2i> 

“ 1^^ l x i“ x 2 ( I x l“ x 2 ^ ^ 


( X>~Y* 2g) *-(iw 


l x l“ x 2 


with 


n ae (i|2, - 


‘a? 3( *l«-*2 a > ( >~V 

^3 r* ? |5 


l x l" x 2 


x l" x 2 


f ' ( l^-Xg | ) 


d 2 f (r ) 
dr 


r “ I x i ” x 2 


f # (|xi-x 2 | } - d 

dr 


r " l x 1 “ x 2 


Defining 


9 g(x--x 2 ) 


(B-18b) 

gives 


(B-19a) 

(B-19b) 
(B-19c ) 

(B-19d) 

(B-20) 
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the dipole-dipole Interaction Eq. (B-14) can be written 


* Va(S <1|2)p a (1) V 2) ' (B - 2l) 

Q£»P 

Eq. (B-19a) is the dipole-dipole interaction law in the 

nonlocal medium. Note that it contains termft in addition 

to the usual free space term, which is obtained by setting 

f equal to one. The extra terms are a consequence of the 

nonlocal nature of the medium. 

Now we want to -examine a special form for e (ic ) 

4 

appropriate to a semiconductor. Walter and Cohen have 
shown that the effects of anisotropy are very small for 
several semiconductors, thus justifying the use of an 
isotropic form for € (it) . For this case, one can define 
a spatial dielectric function €(r) by the relation 


v(r ) 


Ze 

e(r)r * 


(B-22 ) 


where v(r) is the screened Coulomb potential of a point 
charge Ze . In fact, one can show that 


f(r) - , (B-23 ) 

€ (r) 

with f(r) given by Eq. (B-18b). 

5 6 

Using the Penn model, Srinivasan has shown for 
silicon that e(r) essentially reaches its value at 
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infinity by the nearest-neighbor distance. This somewhat 
surprising conclusion can be made plausible as follows. 
With four valence electrons per atom and two atoms per 
primitive cell, the valence electron concentration of the 
diamond structure is 


8 32 

* n V " 75 1 ’ 

o a 

so that the Fermi wavenumber is 


k 


F 


(3ir 2 n y ) 


1/3 



(B-24) 


(B-25a) 


or 


k 


F 


(1.5631853 



For silicon one finds 


(B-25b ) 


-- - 0.5529 A , (B-26 ) 

K F 

which is rather short indeed. One expects screening on 
this length scale . 

To be more quantitative, the results of Walter and 
Cohen 7 can be fit rather closely by the functional form 


e(k) 


1 + 


k 



(B-27) 


where € is the static, bulk dielectric constant. The 
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g 

measured value of € Is 11.7 for silicon. U, ,ing Eq. 
(B-27 ) and performing the Integral In Eq. (B-18b) yields 

f (r) " [7 + (l-|>e-^ . (B-28) 

From the work of Walter and Cohen we find for silicon that 

n - with X - 0.45 , (B-29) 

so that 




pVe - (1.54)^- , 

(B-3 0) 

or a 

screening length of 




1 0 

nVT * °- 56 A • 

(B-31) 

Note 

the close agreement of Eqs. (B-26) and (B-31) . 

Dif- 

ferentiation of Eq . 

(B-28) gives 



f '(r) 

- - Jl u ( 1 * T , 

(B-32a) 


f'(r) 

- eil 2 ( 1 ur . 

(B-32b ) 


In Table Bl f and its first two derivatives are evaluated 
for silicon at the distances of the first few neighbors. 

As can be seen from the Table, f is within one percent 
of its value at infinity by the second neighbor distance. 
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From Eq. (B-19a) we see that the first derivative term has 
the same angular variation as the term with no derivative. 
However, the term with the second derivative has a dif- 
ferent angular variation than the other two terms. This 
could give a sizeable contribution to the dipole-dipole 
interaction energy out to about third nearest neighbors. 
Thus we see that we can write the dipole-dipole interaction 
in the form 


v ii-7 (B ' 33) 

a,P 

+ (exponential terms) , 

where the exponential terms are very short ranged. 


* 


« 
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TABLE B1 

Eqs. (B-28 ) and (B-32) evaluated for the first few 


neighbors 

r 

of silicon. 
f(D 

rf'(r) 

r 2 f '(r) 

r o 

9.932xl0~ 2 

-5 .8 04x1 o” 2 

2 . 432xl0 -1 

r 2 

8.645xl0" 2 

-6. 682x1 0~ 3 

4.572xl0“ 2 

r 3 

8.577xl0“ 2 

-2.405xl0~ 3 

1.93 0x1 0" 2 

r 4 

1.553xl0‘ 2 

-5.554xl0“ 4 

5. 3 74x1 0“ 3 

• 

8.547xl0” 2 

0 

0 
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Appendix C 


QUADRUPOLE MOMENT DUE TO AN ARRAY OF DIPOLES 

We begin with the charge density p(x) due to an 
array of dipoles p ( A k ) located at sites . 

p(x) - - £ J(AK).76(5-it(lK)) . (C-l) 

JtK 

Defining the total quadrupole moment as 

Q ag " J* 3 * x a x 5 p(x) ' (C ’ 2) 

we find after Integration by parts that for the array of 
dipoles 



Now consider the special case of the Raman mode in 

the diamond structure. For this mode u (1,0) - -u (1,1) ■ 

a. ct 

u . Using Eqs. (1.149) and (1.151a) we find that 
a 

P a U,0) - - Sp^ (C-4a) 

P a (i,l) - + Sp^ . (C-4b) 

Clearly there is no total dipole moment due to this mode, 
since the dipole moments in each unit cell are equal and 
opposite and thus add to zero. Using Eq. (C-3 ) and the 
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position vectors Eq. (1.1) we hsve 

<*«* - *»1 Z t [ R a (l>u s + »S <1 > u a] • (c ‘ 5) 

Performing the sum on l and using Eq. (1.5) we obtain 

«„e * 2 »'>i N [ ( ‘ 0 i +6 a2 +# 0 3 >u e + 

+ <‘ei +6 92 +6 9 3 )u «] • (C ' 6) 

where N is the number of unit cells of the crystal. Thus 
we see that the Raman mode has a macroscopic quadrupole 
moment associated with it, where by macroscopic we mean 
proportional to the number of unit cells. Note also that 
only the p^ term of Eq. (1.150) contributes to the quad- 
rupole moment of this mode. 


Appendix D 


THE EWALD METHOD AND ITS APPLICATION 
TO DIPOLE SUMS. 


In this appendix we discuss the Ewald method for 
transforming lattice sums to rapidly convergent form. The 
method is applied to dipole sums in the diamond structure. 
We want to evaluate the sums Eqs. (1.168) 

Q^Oc) - i D'n s U,0|0 > 0)e“ i ^ , ^ U) (D-l ) 

Otp ^ ap 

and 

- ^2 n oS (i,0|0,l)e* ilt,5(i) (D-2) 

A 


with 


fi g UK \l'K') 


6 . 3R UK U'k')r_ UK |x'k') 

a ! E 1 

|?(1K j Z'k') I 3 |K(jJK | X'k ') | 5 

(D-3 ) 


First we define 


A ag <£,x) 


d 


2 




L 

i 


i£*itu) 

e 

[RU)-x 


(D-4) 


Dif ferent iating we obtain 
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, e . c r 


(D-5) 


4 

V s ? !5i [iJr i*r. 


(D-6 ) 


a 2 i 


A a B (k > It) 


f 6 tt a _ 3(R f> (/)-x m ) ( R B U ) -x e ) el g,S (J) 
i |5(i)-x| 3 |3(i)-x| 


(D-7) 


where the prime on the sum indicates omit the l - 0 term. 


Thus we have 


^c4. S ,^ 1 | r ]. r- U.o,-.o,.^*U > 


(D-8 ) 


so that 




n r 

- ^ lim t 
4ir £-0 


A a (-k,x) + 


9x 9x 


x 


(D-9) 


Similarly from Eqs . (D-2) and (D-5) we have 


A „(-£,£<!)) - £ n g(*,0|0,l)e 


-ik-I(X) (D-10) 
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so thut 

- 3 ? A^t-k.Sd)) . (D-ll) 

From Eqs. (D-9) and (D - ■ x > we see that it Is sufficient to 
find an expression for A ^ (ic,x) . Next consider the sum 

_ » _ iit-S(i) 

H(k,x) - 2 2- — CD-12) 

l |J(i)-x| 

using the integral representation 


RU)-x 


J2 

Jtt 


- -4 r do .-P^|R(i)-X 


y d P 
0 


(D-13 ) 


we have 


H(S,x) 



j^(i)-x| 2 e iic*^(X) 


(D-14a) 


or 


H(k,x) 



’£ e “P 2 |R (£) “^| 2 e i ^‘ (&*>-*> 

-l 


(D-14b) 


Def ining 


F(x) .2e'P 2 |^ ) '*| 2 e ii '' ( “ (i) ** ) 
l 


(D-15 ) 
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we note that F(x) - F(x+lt(m)) where lf(m) is a lattice 
translation vector. Thus F(x) is periodic with the 
periodicity of the direct lattice and we can expand it in 
a fourier series* 


F (x) - E F e 1 ^** , (D-16a) 

5 * 

with Fourier coefficients 

F - ~ J* d 3 x FCxJe" 1 ^’* . (D-16b) 

2 u o n 

o 

The sum in Eq. (D-16a) is over all reciprocal lattice 
vectors and the integral in Eq. (E-16b) is over the primi- 
tive unit cell. 

Substituting Eq. (D-15) into Eq. (D-16b) we have 
F, - J-Ej' d 3 x e - l2,i 'e' p2 l 5<1) **| 2 e lS,(5u) - ;) 

g o i n 

^ /n 1 -7 \ 


letting 


y 


X -£(X) 


(D-18 ) 


Eq. (D-17) becomes 



— £ P 

Q 7 •> 

o l w 


.3 

d y 


2 2 

-p y 

e w e 


-i£. 


(R(£)+y ) e «ik*y 


-l 


(D-19 ) 
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where the integral is over the volume of the unit cell at 
the position . Since - 1 we have an 

identical integral for each unit cell so that 

F " 7 r~ £ d 3 y e -P 2 y 2 e -U^)-y (D-20) 

G o V 

where the Integral is now over the volume of the entire 

crystal. We now pass to the limit of a very large crystal 

» 

so that the integral in Eq. (D-20) can be taken over all 
space. Evaluating the integral in Eq. (D-4o) in rectangu- 
lar coordinates we obtain 

, (D _ 21) 

G V 

so that Eq. (D-15) may be written as 


F(x) 


§ e -|^! 2 /4p 2 e i^'x 

n o P G 


(D-22 ) 


Returning to Eq. (D-14b) we break the integral into 
two parts and substitute the Fourier series Eq. (D-22) into 
the first integral to obtain 


H(ic,x) 


^o e 


i<£+£) 


• x 


n — | k+(j | 2 / 4p 

J dp ^ 

0 p 


i 


ik. 


R(i) J dp 

T1 


4 -p 2 |I(f )-X 


(D-23) 


The first integral can be evaluated with the substitution 
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x 



(D-24) 


so that 




-l^lW 

1 e 


(D-25) 


Using the definition of the complimentary error function 


erfc (x) 




(D-26 ) 


the second integral in Eq. (D-23) is 


_2 • -|?(i)-x|j2 

^ e 


erfc(nlit(i).-xl) . 

|!t(i)-x 


(D-27) 


Using Eqs. (D-25) and (D-27) we have from Eq. (D-23) 


H(k,x) 


-l£+5l 2 /4ti 2 


, s erfc(Ti|g(X)-x|) 0 ig-9u) _ ( D -28) 

X |$(A)-x| 

In Eq. (D-28) rapid convergence of both sums can be ob- 
tained by a proper choice of the separation parameter ti . 
Using Eqs. (D-4) and (D-12) we perform the differentiation 
on Eq. (D-28) to obtain 
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( £,J) - fr g .-imlW 

|k+G| 


<*P 




J 2 e lJ ' S( •“ >:'dp [ 2p 2 6 ae -4p 4 <R o ( I ) -x a ) 
* 


-|5u)-x| 2 p 2 


x (R g (i)-Xp)]e 


(D-29) 


using the integral representation 


1 _2 

~ - J dfl e 

i^l v* 


• 2 | •* 1 2 

j ^ .-f I* 

0 


(D-30) 


we find upon differentiation that 


3 2 1 


8 x n dx e X 


2 * 2 4 -|x| 2 p 2 

' 6 oe- 4 P x c x e >e 


(D-31) 


Using Eqs. (D-29) and (D-31) Eq. (D-9) becomes 


&<*> - 


(G -k )(G,-k. ) .* ^,2. J 

O a s L. e l G " k l /4 ti 

li-iTi 2 


n 


+ ^ 375 * 


s - e -lk’R(X)j dp (2p 2 5 oe -4|0 4 R o (i)R e (i)]e'^ <1, l 


o 

+ —m 


*1 

. 2 , 


l- - £5 (^ dp!2p26 ar^ Vp le ' |x|2p2 


dp[2p 2 5 ae -4p 4 x a x p ]e 


. 2^2 
“1*1 P 


(D-32) 


Taking the limit in the last term of Eq. (D-32) we obtain 
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Q «> (S) . .-lSl 2 /4, 2 , ^. ( W< G jrV e -|3-it| 2 /4Ti ! 

k I G 


|5-k 


2 2 


lit • 5( 4 2p 2 6 ae -4p 4 R„( 1 ) R a ( i ) ]e- I R ( 1 > I 0 


+ 7x^7 * 

2 ir l 


nr\ 


a v ' 3 
3 


3ir 


372 o S ’ 


(D-33) 


Where we have explicitly separated out the nonanalytic 
^ > 2 

(j - o term from the sum on G . Note also that the sura 
on L omits the l m 0 term. 

Similarly, using Eqs. (D-ll) and (D-29) we obtain 


(,<?>(£> . e -lS'S(l).-|f|W 

|k| 2 

r'. (G o" k o )(G B" k e ) ) • S<1 )_- |2-k | 2 /4T ! 2 

4- 4 _» _» 3 ^ e © 

3 13-ki 2 

2ir 1 

-|5t(l)-«(l)|V • (D-34) 

x e 


Eqs. (D-33) and (D-34) can be rewritten as 
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M x . k k a i£|2 a 2 (G -k ) (G fl -k a ) i* ^|2.. 2 
. _£Lie'l k l /4t > +Z‘ « « ' | p e~ [ u-k I /4 ti 

iSr “ | 3 - 5 r 


. i y- -ut-fcij-d),, „> "o ’ 1 . 

+ 27? ® : ae u<r|> ^377 s ae <D 35> 


and 


Q ( ?)(S) - e -i*-5(l) -|5| 2 /4n 

aP i ir i " 


:?i 2 . „ 2 


, (0 .-*.>W;-S> 1(3-?). S(X) -|3-£| 2 /4t, 2 

; |3-k 3 


+ 




Tl) 


(D-36) 


where we have defined 


Q 


I ae l)(i ' T ’ ) - «» 


2 » 2 \ 6 


^ 2 2 
.-|S(i)| 0 


(D-3 7a) 


and 


I ^ >U > T ' ) - 77 J 1 dp 

v T| 


a. 


2 ®\r 4 p 4 ( V i)+ V 1)) x 


X <R (*)+R (1)) 


2 2 


-|Jd)+5(i) |p 


(D-3 7b) 


Using the integrals 
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-1 f* t t 2 -aV . 

77 J J dt t • 


1 I , 1 

x -j e““ x +— 3 - erfc(ax) 


77 a' 


2a 


(D-38 ) 


1 m ,22 , „ 3 2 2\ 

4 J dt t V“ 1 - pU +-^—5 e-“ x ) 

fix \2J7 a 4 a I 


+ — r erfc (ax) 

4a 


(D-39 ) 


we can rewrite Eqs . (D-37) as 


Q 


I ifl )(jt » T l ) “ ~f n a g<A,0|0,0)er#c(ri|^(i)|) 


ad 

Q 

o 

+ 7? 


r> A 


crp 


R ( X ) R s ( X ) / n r>\ 

-2— -&-j -l-r^ * 2T1 3 ' 


^JL_ 

|ftu )| 2 |?(i)! z \|ff<x)| 


-T1 2 |R(X)| 2 


(D-40a) 




n 

0 (i,l|0,0)erfc(r 1 |R(2)+^(l) j ) 

* UK 


£0 


*n6 


M. 


(R (X)+R (1))(R (X)+R fl (1)) 
Of a P P 


|it(i)+5(i) 


3n 


|it(X)+R(l) 


2tl‘ 


|ff(X)+£(l) | 2 
-n 2 |£(Z)+£(i) | 2 


(D-40b) 


Eqs. (D-35), (D-36) and (D-40) can now be used to compute 
the sums Eqs. (D-l) and (D-2). All that remains is to 
choose an appropriate value for the separation parameter 



T) . We have found that r\ - gives near optimum con- 

vergence for the diamond structure. 
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Appendix E 


SOME PROPERTIES OF THE DIPOLE CONTRIBUTION 
TO THE DYNAMICAL MATRIX 


In this appendix ve examine some properties of certain 
dipole sums. We begin by considering the dipoles induced 
by the homogeneous deformation 


U a (fK) 


- 2 

e 




d a (K) 


( E-l ) 


where € A are the deformation parameters and d (k) are 

otc a 

the inner displacements. Substituting Eq. (E-l) into 
Eqs. (1.149) and using Eqs. (1.62), (1.4) and (1.5) we 
obtain 


p a (l,0) - £ 2 P^(0,0|6 1 ,l)[d e (l)-d j3 (0)] 

3 i-1 

4 

4 . 2tg 2 p g(0,0|6 1)R (6 1) , (E-2) 

3 Y BY i-1 1 Y i 


and 


p a (f,l) - -p a (i,0) . (E-3) 

Note at this point that th? dipoles induced by the homo- 
geneous deformation are independent of the unit cell index 
so we can write 


p (1,0) - p 

r ot r a 


(E-4) 
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For the important special case of an isotropic homo- 
geneous deformation where 




(E-5) 


and 


d a (K) - 0 , 


(T 




it is easy to verify using Eqs. (1.150) that p^ is ze; 
for this deformation. Thus an isotropic deformation does 
not induce any dipoles and no attendant interaction energy. 
This being the case, it is easy to see why the dipole in- 
teraction does not affect the bulk modulus of the crystal. 

Returning to the more general case and substituting 
Eqs. (E-3 ) and (E-4) into Eqs. (1.153-1.155), we obtain for 
the dipole-dipole interaction energy due to the deformation, 
Eq. (E-l), that 

* dd - z Pci p 3 £;[n ae u,o|i\o)-fi a 3 a,o|r,i)] . 

s & JL 

(E-7) 

In obtaining Eq. (E-7) we have used the fact that 

2^ ag (i,°|i',l) - Zn a& (t,l\l',Q) , (E-8) 

jl jl lit 

which follows from relabeling the sum indices and the pro- 
perty 

^(XkIi'k') - n ag (X'K'|iK) . 
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Defining 


°ae > - J(. n o6 (i,0|l-,0) , CE-10) 

it it 

c^ } - £, n ag (f,o|x',i) , (E-ll) 

we only need show that these terms are zero for an in- 
finite crystal to show that the homogeneous deformation, 
Eq. (E-l . produces no dipole-dipole interaction energy. 
This being the case it then follows that the dipoles do 
not contribute to any of the elastic constants. We now 
proceed to show that and are zero for an in- 

finite crystal. First note that translation invariance 
allows a relabeling of the sums so that 

e'e’ - 2; n aS (ji,o|o,o) <E-i2) 

it ^ 

°<f - £ n a .(i.i|o,o) . (E-13 ) 

Next consider the sum 

L ag (K) -S' (XK |0,0) . (E-14) 

The diagonal terms are easily shown to be zero since using 
Eq. (1.154) 
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L 1X (tt ) 


E' 

i 


fRj(iK)+R^(XK)^(XK)-3Rj(XK) ] 

— ” ' (IukTT 5 


(E-15) 


Using the position vectors Eq. (1.5) and relabeling sum 
indices we note that 

.2, 


E . R?(iK) _ 2 , r|oo _ s . HjjUK > 
t |itUK) I* 1 |S(iK)| 5 1 |3(iK) I 5 


(E-16 ) 


so that we have 


L n (K) - 0 . 


( E-17 ) 


Similar reasoning applies to show that 


L 22 (,t) - Sa* 10 ■ 0 


(E-18) 


Now consider the off-diagonal term 


3R 1 (iK)R 2 (£K) 


L 19 (k) - L - 

i |R( XK ) 


(E-19) 


To show that thi? sum is zero, we can use* a parity argu- 
ment. Noting that we can relabel the summation variables 


by the transformation 

R 1 ( Ak ) - 

R 2 (Xk) - 

R 3 (AK > - 


R 

(E-20a) 

-R . 2 ( A ' K ) 

(E-20b) 

-RgU'tO , 

( E-20c ) 
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we have the new variables 


L* m £ 

*1 *3 


Noting that 


we have that 


i>2 m ** ^ 


l i ■ x i • 


|S(«)| - |S<i'k)| , 


t 12 (K) 


3R. (iK)R (iK) 
2* + t & 

i ifcUK) I 5 


'jo that 


(E-21a) 
(E-21b) 
(E-21c ) 


(E-22) 


(E-23) 


L 12 (K) “ “ L 12 (K) " 0 • (E-24) 

Similar arguments apply to the other off-diagonal elements 
so that 


Since 


L a g(K) - 0 . ( E-20) 

G rf ) "f L ae (0) < E - 21 > 

G IV -2L ae Cl) , 


(E-22) 
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we have that 


<$’ - - ° • <e - 23> 

Thus we find that under the homogeneous deformation, Eq. 
(E-l), 


Since the homogeneous deformation produces no dipole con- 
tribution to the potential energy, one can conclude that 
the elastic constants are not affected by the nonlocal 
dipole-dipole interaction. This is so since a homogeneous 

deformation is sufficient to determine the elastic con- 
2 

stants . 

It is now easy to show that the Raman frequency is 
also unaffected by the nonlocal dipoles. From Eqs. (1.159) 
and (1.161-1.164) we obtain at k - 0 that 

.. 32p 2 

d“(°.°|S*<» -tt?' n a(3 u,ofo,o) 

S A 

32p i - 

- l-jJ-S n a3 (i,l|0,0) , (E-24) 

S i 

and 


d“(o,o|E- 


°) - -d“(0,i|S-0) 


(E-25 ) 


Using Eqs. (E-14) and (E-20) we have 
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D ^(KK'|S-°) - 0 . (E-26) 

Thus we see that the dipoles do not contribute to the Raman 
frequency. 

Similar types of symmetry arguments ran be used to 
show the properties stated in Chapter 1, though the details 
are too length> .o present here. 
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Appendix F 


PROOF THAT S (1) (iK) AND * (2 ) Uk) ARE ZERO 
a a 

IN THE DIAMOND STRUCTURE 

We begin by showing that 4^*^ (IK) has the same trans- 

a 

formation properties as the first order atomic force con- 
stants. This is obviously true since an isotropic homogen- 
eous deformation does not change the symmetry of the cry- 
stal . 

Using the definition Eq. (2.23c) 

S^ 1 ) (LK) - Z * (LK|l/K')R-(l/K') , (F-l ) 


we apply the transformation law Eq. (1.35) for the second 
order force constants under the symmetry operation Eq, 
(1.34) to obtain 


8 (1) (LX) 
a 


L tc 3 


Z: S fl , 4 , AM \A'x') 

^ au \iU 


Z S pc R c (i'tc')+v p (S )+Rg (m) 


(F-2) 


or 


4 (1) (LK) - Z Z S S Q ,S 0 * tJ (iK |Z'K')R U\') 
a Ji'K'S aUn &V to 1 a 


l K 0 y,l/c 


+ 1 - K'e £ s ^ s 3 ‘'V'' Uk| ‘' < ' ,[v 9 <e,+r 9 ( ' b)I 


(F-3 ) 
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Using the Infinitesimal translation Invariance condition 


1 


£ l_ a UK| i\') - 0 , 


i'k' 


ere 


(r-4) 


Eq . (F-3 ) becomes 


»^ 1) (LK) 


£ £s S 0i S fl * ,(iK \l\ # )R d'K') . 
i'K'3 Wa 3l/ 00 ^ a 


(F-5) 


Now we note that since S is a real orthogonal matrix 


£ S 

e 


e^ s eo 


£ s 
e 


" 1 c 



(F-6 ) 


therefore 


4 (1) (LK) 

a 


£ 

I'k' 


£ S i, ! , (iK |X'k')R U'k') , 

<*» VO V.V 1 a 

(F-7 ) 


or 


4 (1) (LK) - £ S 

Ct C 


£ 4 (IkU'k^r U'O 

^ l\'v » v v 


(F-8) 


By the definition Eq. (F-l) 


- 2 4 UK|i'K')R U'0 , 

^ X'K'l/ ^ V 


(F-9) 


so that 
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* ( 1 ) (LK)-Es t ( 1 ) (iK) . (F-10) 

at au |i 

Eq. (F-10) is the same transformation law obeyed by the 

2 

first order atomic force constants. Similar reasoning 

(2) 

can be used to show that '(Ik) as defined by Eq. 
(2.23d) transforms as 

* (2) (LK) - E S i (2) (lK) . (F-ll) 

a oty, c 

Now that we have established the transformation law, 
we apply the inversion operation Eqs. (1.36) and (1.37) 
to obtain 

»^ 1) (-1,1) - - * ( ^ 1) U, 0) . (F-12 ) 

Lattice translation invariance gives 

»^ 1) U+m,K) - i^hl.K) , (F-13) 

so that the coefficients are independent of the unit cell 
index 


* (1) (iK) - t (1) (0,K) . (F-14) 

ot at 

combining Eqs. (F-12) and (F-14) we have 

*^ 1) (1,0) - *^ 1) (0,0) - - 4 ( ^ 1) (0,1) - - i^U,!) . 


(F-15 ) 
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The simple two-fold rotation 


element 


3 




gives 




^ 1 } ( 0 , 0 ) - 0 , 


and 


*i 1) (0,0) - 


- »^ X) (0,0) - 0 


The further operation 6 


2y 


2y 



gives 


( 0 , 0 ) 


- *^ 1 ) ( 0 , 0 ) - 0 


(F-16) 


(F-17a) 


(F-17b ) 


(F-18 ) 


(F-19 ) 


Eqs. (F-17) and (F-19) together with Eq. (F-15) show that 


* (1> (Xk) - 0 . 

Of 


(F * )) 



(2 ) 

Since (4k) obeys the same transformation law, 

it follows that it also is zero in .he diamond structure. 
Note that the proof only applies to an infinite crystal, 
since surfaces break the translation invariance condition 
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